# Optimal (Control of) Intervention Strategies for Malaria Epidemic in Karonga District, Malawi

Authors: Peter M. Mwamtobe; Shirley Abelman; J. Michel Tchuenche; Ansley Kasambara
Journal: Abstract and Applied Analysis
 (2014)

---

## Abstract

Malaria is a public health problem for more than 2 billion people globally. About 219
million cases of malaria occur worldwide and 660,000 people die, mostly (91%) in the African
Region despite decades of efforts to control the disease. Although the disease is preventable, it
is life-threatening and parasitically transmitted by the bite of the femaleAnopheles mosquito.
A deterministic mathematical model with intervention strategies is developed in order to
investigate the effectiveness and optimal control strategies of indoor residual spraying (IRS),
insecticide treated nets (ITNs) and treatment on the transmission dynamics of malaria in
Karonga District, Malawi. The effective reproduction number is analytically computed,
and the existence and stability conditions of the equilibria are explored. The model does
not exhibit backward bifurcation. Pontryagin’s Maximum Principle which uses both the
Lagrangian and Hamiltonian principles with respect to a time dependent constant is used to
derive the necessary conditions for the optimal control of the disease. Numerical simulations
indicate that the prevention strategies lead to the reduction of both the mosquito population
and infected human individuals. Effective treatment consolidates the prevention strategies.
Thus, malaria can be eradicated in Karonga District by concurrently applying vector control via ITNs and IRS complemented with timely treatment of infected
people.


---

## Body

## 1. Introduction

Malaria is a vector-borne infectious disease found mainly in tropical regions (Sub-Saharan Africa, Central and South America, the Indian subcontinent, Southeast Asia, and the Pacific islands) [1]. It is a life-threatening disease transmitted through the bites of infected mosquitoes [2]. There are four different types ofPlasmodium parasites:Plasmodium falciparum (the only parasite which causes malignant malaria),Plasmodium vivax (causes benign malaria with less severe symptoms; the vector can remain in the liver for up to three years and can lead to a relapse),Plasmodium malariae (also causes benign malaria and is relatively rare), andPlasmodium ovale (causes benign malaria and can remain in the blood and liver for many years without causing symptoms). This study focuses mainly on malignant malaria. Severe malaria can affect the patient’s brain and central nervous system and can be fatal [2, 3]. Furthermore, Medicinenet [4] has reported another relatively new speciesPlasmodium knowlesi which has been causing malaria in Malaysia and areas of Southeast Asia. It is also a dangerous species that is typically found only in long-tailed and pigtail macaque monkeys. LikeP. falciparum,P. knowlesi may be deadly to anyone infected [5].Beyond the human toll, malaria wreaks significant economic havoc in endemic regions, decreasing gross domestic product (GDP) by as much as 1.3% in countries with high levels of transmission—the disease accounts for up to 40% of public health expenditures, 30–50% of in-patient hospital admissions, and up to 60% of out-patient health clinic visits [2]. The Government of Malawi has put in place several control strategies through the National Malaria Control Programme to reduce and possibly (ultimately) eliminate malaria. The main strategic areas that have been identified for scaling-up of malaria control activities include malaria case management, intermittent preventive treatment (IPT) of pregnant women using sulfadoxine-pyrimethamine (SP), and malaria prevention with special emphasis on the use of ITNs as well as IRS [8].Various compartmental models for the spread of malaria have been proposed [10–12]. Extensive review of mathematical models of malaria dynamics using S
E
I
R
S-type models and their variants can be found in Chiyaka et al. [13, 14]. For information on the unexpected stability of malaria elimination and eradication and additional literature on malaria dynamics see Smith et al. [15], who also showed that if mosquito birth rate is increased by 25%, then the minimal effective spraying period is reduced by half, and if doubled, the period is reduced by three-quarters. Our goal is to assess the role of optimal control of preventive measures, namely, ITNs and IRS as well as treatment of the transmission dynamics of malaria in Karonga District, Malawi, without actually targeting a certain high-risk group. Despite a plethora of studies on the dynamics of malaria and its control (see Chiyaka et al. [14] and Smith et. al. [15]), to the best of our knowledge, the proposed model with exposed immigrants is seemingly new: the use of optimal values of a combination of three intervention strategies, namely, indoor residual spraying (IRS), insecticide treated nets (ITNs), and treatment to investigate the effectiveness and optimal control strategies in the transmission dynamics of malaria in a locality in a resource constrained setting.


## 2. Model Formulation and Analysis

We formulate an optimal control model for malaria with the population under study being subdivided into compartments according to individuals’ disease status. We consider the total population sizes denoted byN
h
(
t
) and N
v
(
t
) for the human hosts andAnopheles female mosquitoes, respectively. We employ the S
E
I
R
S framework to describe a disease with temporary immunity on recovery from infection. The S
E
I
R
S model indicates that the passage of individuals is from the susceptible class, S
h, to the exposed class, E
h, then to the infectious class, I
h, and finally to the recovery class, R
h. S
h
(
t
) represents the number of individuals not yet infected with the malaria parasite at time t. The latent or exposed class E
h
(
t
) represents individuals who are infected but not yet infectious. Individuals in the I
h
(
t
) class are infected with malaria and are capable of transmitting the disease to susceptible mosquitoes. R
h
(
t
) represents the class of individuals who have temporarily recovered from the disease. The susceptible human population is increased by recruitment (birth) at a constant rate, Λ
h, while others are generated through migration by (
1
-
κ
1
)
θ
N
h, where κ
1 is the proportion of infected immigrants into the exposed class and θ is the rate at which people migrate into Karonga District. All the recruited individuals are assumed to be naive when joining the community.There is some finite probability,β
v
h, that the parasites (in the form of sporozoites) will be passed onto the humans when an infectious femaleAnopheles mosquito bites a susceptible human. The parasite then moves to the liver where it develops into its next life stage, merozoites. Using the approach adopted in [6], susceptible individuals acquire malaria through contact with infectious mosquitoes at the rate ϑ. The infected person moves to the exposed class at the rate (
1
-
u
1
)
λ
h
I
v
S
h. The preventive variable u
1
(
t
)
∈
[
0
,
1
] represents the use of ITNs as a means of minimizing or eliminating mosquito-human contacts. After a certain period of time, the parasite (in the form of merozoites) enters the blood stream, usually signaling the clinical onset of malaria. Then the exposed individuals become infectious and progress to the infected state at a constant rate α
h. The individuals who have experienced infection may recover with temporary immunity at a constant rate ρ and move to the recovery class, while some infectious humans after recovery without immunity become immediately susceptible again at the rate (
1
-
ρ
). Infectious individuals recover due to treatment at a rate η
u
2 with u
2
(
t
)
∈
[
0
,
1
] representing the control effort on treatment and ϕ being the proportion of individuals who recover spontaneously. Recovered individuals lose immunity at a rate ψ. The natural and disease induced death rates are μ
h and δ
h, respectively. The disease induced death rate is very small in comparison with the recovery rate.The mosquito populationN
v is divided into three compartments: susceptible S
v
(
t
); exposed E
v
(
t
); and infectious I
v
(
t
). FemaleAnopheles mosquitoes enter the susceptible class through birth at a rate Λ
v. The parasites in the form of gametocytes enter the mosquito population with probability β
h
v. This happens when the mosquito bites an infectious human and the mosquito moves from the susceptible to the exposed class. Mosquitoes are assumed to suffer death due to natural causes at a rate μ
v. The exposed mosquitoes progress to the class of symptomatic mosquitoes I
v at a rate α
v. It is assumed that the disease does not induce death to the mosquito population. Finally, the mortality rate of the mosquito population increases at a rate proportional to τ
u
3
(
t
), where τ
>
0 is a rate constant when using IRS and u
3
(
t
)
∈
[
0
,
1
] is a constant rate of IRS. Figure 1 represents a schematic flow diagram of the proposed model.Figure 1
The malaria model with interventions flowchart.The model state variables are represented in Table1. Table 2 represents prevention and control strategies practised in the district. Table 3 shows parameters of the model.Table 1
State variables of the malaria model.
Symbol
Description
S
h
(
t
)
Number of susceptible individuals at timet
E
h
(
t
)
Number of exposed individuals at timet
I
h
(
t
)
Number of infectious humans at timet
R
h
(
t
)
Number of recovered humans at timet
S
v
(
t
)
Number of susceptible mosquitoes at timet
E
v
(
t
)
Number of infected mosquitoes at timet
I
v
(
t
)
Number of infectious mosquitoes at timet
N
h
(
t
)
Total number of individuals at timet
N
v
(
t
)
Total mosquito population at timetTable 2
Prevention and control variables in the model.
Symbol
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 spraying (IRS)
τ
Rate constant due to use of indoor residual spraying
η
Rate constant due to use of treatment effortTable 3
Parameters variables of the malaria model.
Symbol
Description
Λ
h
Recruitment rate of individuals by birth
κ
1
Proportion of exposed immigrants into exposed class
θ
Proportion of people migrating into Karonga District
μ
h
Per capita natural death rate of humans
μ
v
Per capita natural death rate of mosquitoes
δ
h
Per capita disease induced death rate for humans
α
h
Progression rate of humans from exposed state to the infectious state
β
v
h
Probability that a bite results in transmission of infection to human
β
h
v
Probability that a bite results in transmission of the parasite from an infectious human to the susceptible mosquitoes
ϑ
Biting rate of mosquito
λ
h
Force of infection for susceptible humans to exposed individuals
λ
v
Force of infection for susceptible mosquitoes to exposed mosquito class
ρ
Per capita rate of recovery with temporary immunity
ϕ
Proportion of spontaneous individual recovery
ψ
Per capita rate of loss immunity
α
v
Progression of exposed mosquitoes into infected mosquitoesThe state variables in Table1, the prevention and control parameters in Table 2, and the model parameters in Table 3 for the malaria model satisfy (1). It is assumed that all state variables and parameters of the model which monitors human and mosquito populations are positive for all t
≥
0. We will, therefore, analyse the model in a suitable region.The assumptions lead to the following deterministic system of nonlinear ordinary differential equations which describe the evolutionary dynamics of a malaria model with a combination of interventions:(1)
d
S
h
d
t
=
Λ
h
+
(
1
-
κ
1
)
θ
N
h
+
(
ϕ
+
η
u
2
)
(
1
-
ρ
)
I
h
-
(
1
-
u
1
)
λ
h
S
h
-
μ
h
S
h
+
ψ
R
h
,
d
E
h
d
t
=
(
1
-
u
1
)
λ
h
S
h
+
κ
1
θ
N
h
-
α
h
E
h
-
μ
h
E
h
,
d
I
h
d
t
=
α
h
E
h
-
(
ϕ
+
η
u
2
)
(
1
-
ρ
)
I
h
-
(
ϕ
+
η
u
2
)
ρ
I
h
-
(
μ
h
+
δ
h
)
I
h
,
d
R
h
d
t
=
(
ϕ
+
η
u
2
)
ρ
I
h
-
(
μ
h
+
ψ
)
R
h
,
d
S
v
d
t
=
Λ
v
-
λ
v
S
v
-
(
μ
v
+
τ
u
3
)
S
v
,
d
E
v
d
t
=
λ
v
S
v
-
α
v
E
v
-
(
μ
v
+
τ
u
3
)
E
v
,
d
I
v
d
t
=
α
v
E
v
-
(
μ
v
+
τ
u
3
)
I
v
,
where λ
h
=
β
v
h
ϑ
I
v
/
N
h, λ
V
=
β
h
v
ϑ
I
h
/
N
h.The termβ
v
h
ϑ
S
h
I
v
/
N
h denotes the rate at which the human hosts S
h become infected by infectious mosquitoes I
v and β
h
v
ϑ
S
v
I
h
/
N
h refers to the rate at which the susceptible mosquitoes S
v are infected by the infectious human hosts I
h.Adding the first six equations of model (1) and assuming that there is no disease induced death, that is, δ
h
=
0, gives d
N
h
/
d
t
=
Λ
h
-
μ
h
N
h, so that N
h
(
t
)
→
Λ
h
/
μ
h as t
→
∞ [16]. Thus, Λ
h
/
μ
h is an upper bound of N
h
(
t
) provided that N
h
(
0
)
≤
Λ
/
μ
h. Further, if N
h
(
0
)
>
Λ
h
/
μ
h, then N
h
(
t
) will decrease to this level, Λ
h
/
μ
h. Similar calculation for the vector equations shows that N
v
→
Λ
v
/
μ
v as t
→
∞.The feasible region(2)
Φ
=
{
Λ
v
μ
v
(
S
h
,
E
h
,
I
h
,
R
h
,
S
v
,
E
v
,
I
v
)
∈
R
+
7
:
N
h
≤
Λ
h
μ
h
=
N
h
*
,
N
v
≤
Λ
v
μ
v
=
N
v
*
}
is positive-invariant and attracting.Lemma 1.
The regionΦ
∈
R
+
7 is positively invariant for the model system (1) with initial conditions in R
+
7.Proof.
Lett
¯
=
sup
⁡
{
t
>
0
:
S
h
>
0
,
E
h
>
0
,
I
h
>
0
,
R
h
>
0
,
S
v
>
0
,
E
v
>
0
,
I
v
>
0
}
∈
[
0
,
t
] give t
¯
>
0. The first equation of model (1) gives
(3)
d
S
h
d
t
=
Λ
h
+
(
1
-
κ
1
)
θ
N
h
+
(
ϕ
+
η
u
2
)
(
1
-
ρ
)
I
h
-
(
1
-
u
1
)
λ
h
S
h
-
μ
h
S
h
+
ψ
R
h
≥
Λ
h
-
(
(
1
-
u
1
)
λ
h
+
μ
h
)
S
h
,
which can be rewritten as
(4)
d
d
t
[
S
h
(
t
)
e
∫
0
t
‍
(
1
-
u
1
)
λ
h
(
s
)
d
s
+
μ
h
t
]
≥
Λ
h
e
∫
0
t
‍
(
1
-
u
1
)
λ
h
(
s
)
d
s
+
μ
h
t
.
Therefore,
(5)
S
h
(
t
¯
)
e
∫
0
t
¯
‍
{
(
1
-
u
1
)
λ
h
(
s
)
d
s
}
+
μ
h
t
¯
-
S
h
(
0
)
≥
∫
0
t
¯
Λ
h
e
∫
0
u
‍
(
1
-
u
1
)
λ
h
(
w
)
d
w
+
u
h
u
d
u
,
so that
(6)
S
h
(
t
¯
)
≥
S
h
(
0
)
e
-
(
∫
0
t
¯
‍
(
1
-
u
1
)
λ
h
(
s
)
d
s
+
μ
h
t
¯
)
+
e
-
(
∫
0
t
¯
‍
(
1
-
u
1
)
λ
h
(
s
)
d
s
+
μ
h
t
¯
)
×
{
∫
0
t
¯
‍
Λ
h
e
∫
0
u
‍
(
1
-
u
1
)
λ
h
(
w
)
d
w
+
μ
h
u
d
u
}
>
0
.
Similarly, it can be shown that E
h
>
0, I
h
>
0, R
h
>
0, S
v
>
0, E
v
>
0, and I
v
>
0, for all t
>
0. This completes the proof.


### 2.1. Existence and Stability of Equilibrium Points

We analyse system (1) to obtain the equilibrium points of the system and their stability. Let (
S
h
*
,
E
h
*
,
I
h
*
,
R
h
*
,
S
v
*
,
E
v
*
,
I
v
*
) be the equilibrium points of system (1). At an equilibrium point, we have
(7)
S
h
′
(
t
)
=
E
h
′
(
t
)
=
I
h
′
(
t
)
=
R
h
′
(
t
)
=
S
v
′
(
t
)
=
E
v
′
(
t
)
=
I
v
′
(
t
)
=
0
.


#### 2.1.1. Disease-Free Equilibrium,E
0

In the absence of malaria, that is,E
h
*
=
I
h
*
=
R
h
*
=
E
v
*
=
I
v
*
=
0, model system (1) has an equilibrium point called the disease-free equilibrium, E
0, and is given by
(8)
E
0
=
(
Λ
h
μ
h
,
0,0
,
0
,
Λ
h
μ
v
,
0,0
)
.
To establish the linear stability of E
0, we employ van den Driessche and Watmough’s next generation matrix approach [17]. A reproduction number obtained this way determines the local stability of the disease-free equilibrium point for R
e
<
1 and instability for R
e
>
1. Following van den Driessche and Watmough [17], the associated next generation matrices F and V of system (1) can be determined from F
i and V
i, respectively, where
(9)
F
i
=
[
(
1
-
u
1
)
β
v
h
ϑ
I
v
S
h
N
h
+
κ
1
θ
N
h
0
β
h
v
ϑ
I
h
S
v
N
h
0
]
,
V
i
=
[
(
α
h
+
μ
h
)
E
h
(
ϕ
+
η
u
2
+
μ
h
+
δ
h
)
I
h
-
α
h
E
h
(
α
v
+
μ
v
+
τ
u
3
)
E
v
(
μ
v
+
τ
u
3
)
I
v
-
α
v
E
v
]
.
The Jacobian matrix of F
i and V
i is given by
(10)
F
=
[
0
0
0
r
1
0
0
0
0
0
j
0
0
0
0
0
0
]
,
V
=
[
m
0
0
0
-
α
h
n
0
0
0
0
a
0
0
0
-
α
v
b
]
,
where
(11)
r
1
=
(
1
-
u
1
)
β
v
h
ϑ
,
j
=
β
h
v
ϑ
Λ
v
μ
h
Λ
h
μ
v
,
m
=
α
h
+
μ
h
,
n
=
ϕ
+
η
u
2
+
μ
h
+
δ
h
,
a
=
α
v
+
μ
v
+
τ
u
3
,
b
=
μ
v
+
τ
u
3
.
Algebraic manipulation of the matrices leads to the effective reproduction number
(12)
R
e
=
(
×
(
α
v
+
μ
v
+
τ
u
3
)
(
μ
v
+
τ
u
3
)
μ
v
Λ
h
ϑ
2
)
-
1
β
h
v
β
v
h
ϑ
2
α
h
α
v
(
1
-
u
1
)
Λ
v
μ
h
f
f
i
×
(
ϑ
2
(
α
h
+
μ
h
)
(
ϕ
+
η
u
2
+
μ
h
+
δ
h
)
f
f
f
i
f
i
×
(
α
v
+
μ
v
+
τ
u
3
)
(
μ
v
+
τ
u
3
)
μ
v
Λ
h
ϑ
2
)
-
1
)
1
/
2
,
where
(13)
R
e
v
=
β
h
v
α
v
ϑ
Λ
v
μ
v
(
α
v
+
μ
v
+
τ
u
3
)
(
μ
v
+
τ
u
3
)
is the contribution of the mosquito population when it infects the humans, and
(14)
R
e
h
=
β
v
h
ϑ
μ
h
α
h
(
1
-
u
1
)
Λ
h
(
α
h
+
μ
h
)
(
ϕ
+
η
u
2
+
μ
h
+
δ
h
)
is the human contribution when they infect the mosquitoes.The expression for the effective reproduction number,R
e, has a biological meaning that is readily interpreted from terms under the square root sign. Consider the following terms.(i)
β
v
h
ϑ
α
v
/
(
α
v
+
μ
v
+
τ
u
3
)
(
μ
v
+
τ
u
3
) represents the number of secondary human infections caused by one infected mosquito vector.
(ii)
β
h
v
ϑ
α
h
/
(
α
h
+
μ
h
)
(
ϕ
+
η
u
2
+
μ
h
+
δ
h
) represents the number of secondary mosquito infections caused by one infected human host.The square root represents the geometric mean of the average number of secondary host infections produced by one vector and the average number of secondary vector infections produced by one host. This effective reproduction number serves as an invasion threshold both for predicting outbreaks and evaluating control strategies that would reduce the spread of the disease. The threshold quantity,R
e, measures the average number of secondary cases generated by a single infected individual in a susceptible human population [18], where a fraction of the susceptible human population is under prevention and the infected class is under treatment. In the absence of any protective measure, the effective reproduction number R
e with treatment is
(15)
R
e
t
=
β
h
v
β
v
h
ϑ
2
α
h
α
v
Λ
v
μ
h
(
α
h
+
μ
h
)
(
ϕ
+
η
u
2
+
μ
h
+
δ
h
)
(
α
v
+
μ
v
)
μ
v
2
Λ
h
.
Also if ITNs are the only intervention strategy, then
(16)
R
e
N
=
β
h
v
β
v
h
ϑ
2
α
h
α
v
Λ
v
μ
h
(
1
-
u
1
)
(
α
h
+
μ
h
)
(
ϕ
+
μ
h
+
δ
h
)
(
α
v
+
μ
v
)
μ
v
2
Λ
h
.
Similarly, if IRS is the only means of protection, then
(17)
R
e
s
=
(
β
h
v
β
v
h
ϑ
2
α
h
α
v
μ
h
Λ
v
fffif
×
(
(
α
h
+
μ
h
)
(
ϕ
+
μ
h
+
δ
h
)
(
α
v
+
μ
v
+
τ
u
3
)
f
f
f
f
f
f
×
(
μ
v
+
τ
u
3
)
μ
v
Λ
h
)
-
1
)
1
/
2
.
The value of the basic reproduction number can be obtained from the value of the effective reproduction number when there are no control measures (
u
1
=
u
2
=
u
3
=
0
):
(18)
R
0
=
β
h
v
β
u
h
ϑ
2
α
h
α
v
Λ
v
μ
h
(
α
h
+
μ
h
)
(
ϕ
+
μ
h
+
δ
h
)
(
α
v
+
μ
v
)
μ
v
2
Λ
h
.
In general, it is easy to prove that
(19)
R
e
≤
R
0
for 0
≤
u
1, u
3
≤
1, due to reduction of likelihood of infection by protection. This implies that ITNs and IRS have a positive impact on the malaria dynamics as they contribute to the reduction of secondary infections. Therefore, from van den Driessche and Watmough [17] (Theorem 2), the following result holds.Lemma 2.
The disease-free equilibrium,E
0, of the malaria model with intervention strategies (1), given by (8), is locally asymptotically stable if R
e
<
1 and unstable if R
e
>
1.


#### 2.1.2. Endemic Equilibrium Point,E
1

When malaria is present, the model system (1) has a steady state, E
1, called the endemic equilibrium. In order to establish the stability of E
1, we express system (1) in dimensionless variables, namely, S
h
/
N
h
=
s, E
h
/
N
h
=
e, I
h
/
N
h
=
i, R
h
/
N
h
=
r, S
v
/
N
v
=
x, E
v
/
N
v
=
y, and I
v
/
N
v
=
z, where d
N
h
/
d
t
=
Λ
h
-
μ
h
N
h
-
δ
h
I
h, d
N
v
/
d
t
=
Λ
v
-
μ
v
N
v. Then differentiating with respect to time, t, respectively results in the following system:
(20)
d
s
d
t
=
1
N
h
[
d
S
h
d
t
-
s
d
N
h
d
t
]
=
1
N
h
[
Λ
h
+
(
1
-
κ
1
)
θ
N
h
+
(
ϕ
+
η
u
2
)
(
1
-
ρ
)
i
N
h
f
f
f
f
i
f
-
(
1
-
u
1
)
β
v
h
ϑ
i
s
N
h
-
μ
h
s
N
h
+
ψ
r
N
h
]
.
Let m
1
=
(
1
-
κ
1
)
θ, n
1
=
(
ϕ
+
η
u
2
)
(
1
-
ρ
), and q
=
(
1
-
u
1
)
β
v
h
ϑ.Then(21)
d
s
d
t
=
1
N
h
[
Λ
h
+
m
1
N
h
+
n
1
i
N
h
-
q
i
s
N
h
-
u
h
s
N
h
+
ψ
r
N
h
]
-
s
N
h
[
Λ
h
-
δ
h
i
N
h
-
(
μ
h
-
θ
)
N
h
]
=
Λ
h
N
h
-
[
Λ
h
N
h
+
θ
-
δ
h
i
]
s
+
(
1
-
κ
1
)
θ
+
(
ϕ
+
η
u
2
)
(
1
-
ρ
)
i
-
(
1
-
u
1
)
β
v
h
ϑ
z
s
+
ψ
r
,
d
e
d
t
=
1
N
h
[
d
E
d
t
-
e
d
N
h
d
t
]
=
1
N
h
[
(
1
-
u
1
)
β
v
h
ϑ
z
s
N
h
+
κ
1
θ
N
h
-
α
h
e
N
h
-
μ
h
e
N
h
]
-
e
N
h
[
Λ
h
-
δ
h
i
N
h
-
(
μ
h
-
θ
)
N
h
]
=
(
1
-
u
1
)
β
v
h
ϑ
z
s
+
κ
1
θ
-
[
Λ
h
N
h
-
δ
h
i
+
α
h
+
θ
]
e
,
d
i
d
t
=
1
N
h
[
d
I
d
t
-
i
d
N
h
d
t
]
=
1
N
h
[
α
h
e
N
h
-
(
ϕ
+
η
u
2
)
(
1
-
ρ
)
i
N
h
fffffffffi
-
(
ϕ
+
η
u
2
)
ρ
i
N
h
-
(
μ
h
+
δ
h
)
i
N
h
]
-
i
N
h
[
Λ
h
-
δ
h
i
N
h
-
(
μ
h
-
θ
)
N
h
]
=
α
h
e
-
[
Λ
h
N
h
+
ϕ
+
η
u
2
+
δ
h
+
θ
-
δ
h
i
]
i
,
d
r
d
t
=
1
N
h
[
d
R
d
t
-
r
d
N
h
d
t
]
=
1
N
h
[
(
ϕ
+
η
u
2
)
ρ
i
N
h
-
(
μ
h
+
ψ
)
r
N
h
]
-
r
N
h
[
Λ
h
-
δ
h
i
N
h
-
(
μ
h
-
θ
)
N
h
]
=
(
ϕ
+
η
u
2
)
ρ
i
-
[
Λ
h
N
h
+
ψ
+
θ
-
δ
h
i
]
r
,
d
x
d
t
=
1
N
v
[
d
S
v
d
t
-
x
d
N
v
d
t
]
=
1
N
v
[
Λ
v
-
β
h
v
ϑ
i
x
N
v
-
(
μ
v
+
τ
u
3
)
x
N
v
]
-
x
N
v
[
Λ
v
-
μ
v
N
v
]
=
[
Λ
v
N
v
-
x
Λ
v
N
v
]
-
β
h
v
ϑ
i
x
-
τ
u
3
x
,
d
y
d
t
=
1
N
v
[
d
E
v
d
t
-
y
d
N
v
d
t
]
=
1
N
v
[
β
h
v
ϑ
i
x
N
v
-
α
v
y
N
v
-
(
μ
v
+
τ
u
3
)
y
N
v
]
-
y
N
v
[
Λ
v
-
μ
v
N
v
]
=
β
h
v
ϑ
i
x
-
[
Λ
v
N
v
+
α
v
+
τ
u
3
]
y
,
d
z
d
t
=
1
N
v
[
d
I
v
d
t
-
Z
d
N
v
d
t
]
=
1
N
v
[
α
v
y
N
v
-
(
μ
v
+
τ
u
3
)
z
N
v
]
-
z
N
v
[
Λ
v
-
μ
v
N
v
]
=
α
v
y
-
[
Λ
v
N
v
+
τ
u
3
]
z
,
d
N
h
d
t
=
[
Λ
h
N
h
-
μ
h
-
δ
h
i
]
N
h
,
d
N
v
d
t
=
[
Λ
v
N
v
-
μ
v
]
N
v
.
The system can now be reduced to a nine-dimensional system by eliminating s and x since s
=
1
-
e
-
i
-
r and x
=
1
-
y
-
z. The following feasible region
(22)
Φ
1
=
{
Λ
v
μ
v
e
,
i
,
r
,
N
h
,
y
,
z
,
N
v
∈
R
+
7
∣
e
≥
0
,
i
≥
0
,
r
≥
0
,
e
+
i
+
r
≤
1
,
N
h
≤
Λ
h
μ
h
,
y
≥
0
,
z
≥
0
,
y
+
z
≤
1
,
N
v
≤
Λ
v
μ
v
}
can be shown to be positively invariant. R
+
7 denotes the nonnegative cone of R
7 including its lower dimensional faces. Thus we have the following system of equations:
(23)
d
e
d
t
=
(
1
-
u
1
)
(
1
-
e
-
i
-
r
)
β
v
h
ϑ
z
+
κ
1
θ
-
[
Λ
h
N
h
-
δ
h
i
+
α
h
+
θ
]
e
,
d
i
d
t
=
α
h
e
-
[
Λ
h
N
h
+
ϕ
+
η
u
2
+
δ
h
+
θ
-
δ
h
i
]
i
,
d
r
d
t
=
(
ϕ
+
η
u
2
)
ρ
i
-
[
Λ
h
N
h
+
ψ
+
θ
+
-
δ
h
i
]
r
,
d
N
h
d
t
=
Λ
h
-
δ
h
i
N
h
-
μ
h
N
h
,
d
y
d
t
=
(
1
-
y
-
z
)
β
h
v
ϑ
i
-
[
Λ
v
N
v
+
α
v
+
τ
u
3
]
y
,
d
z
d
t
=
α
v
-
[
Λ
v
N
v
+
τ
u
3
]
z
,
d
N
v
d
t
=
Λ
v
-
μ
v
N
v
.
We seek to establish whether a unique endemic equilibrium exists. This is done by making more realistic assumptions, namely, that the protective control measures may not be totally effective.Existence and Uniqueness of Endemic Equilibrium,
E
1. To compute the steady states of the system (23), we set the derivative with respect to time in (23) equal to zero, and after simplification the following algebraic equations are obtained:
(24)
(
1
-
u
1
)
(
1
-
e
-
i
-
r
)
β
v
h
ϑ
z
+
κ
1
θ
=
[
Λ
h
N
h
-
δ
h
i
+
α
h
+
θ
]
e
,
α
h
e
=
[
Λ
h
N
h
+
ϕ
+
η
u
2
+
δ
h
+
θ
-
δ
h
i
]
i
,
(
ϕ
+
η
u
2
)
ρ
i
=
[
Λ
h
N
h
+
ψ
+
θ
-
δ
h
i
]
r
,
Λ
h
N
h
=
δ
h
i
+
μ
h
,
(
1
-
y
-
z
)
β
h
v
ϑ
i
=
[
Λ
v
N
v
+
α
v
+
τ
u
3
]
y
,
α
v
y
=
[
Λ
v
N
v
+
η
u
3
]
z
,
Λ
v
N
v
=
μ
v
.
To calculate the dimensionless proportions in terms of i, consider the first equation in the system (24):
(25)
(
1
-
i
-
r
)
p
z
-
p
z
e
+
κ
1
θ
=
[
μ
h
+
α
h
]
e
,
(
μ
h
+
α
h
)
e
+
(
1
-
u
1
)
β
v
h
ϑ
z
e
=
(
1
-
u
1
)
(
1
-
i
-
r
)
β
v
h
ϑ
z
+
κ
1
θ
,
(
μ
h
+
α
h
+
(
1
-
u
1
)
β
v
h
z
)
e
=
(
1
-
u
1
)
(
1
-
i
-
r
)
β
v
h
ϑ
z
+
κ
1
θ
,
e
=
(
1
-
u
1
)
(
1
-
i
-
r
)
β
v
h
ϑ
z
+
κ
1
θ
μ
h
+
α
h
+
(
1
-
u
1
)
β
v
h
ϑ
z
,
where p
=
(
1
-
u
1
)
β
v
h
ϑ.From the second equation we have(26)
α
h
e
=
(
μ
h
+
ϕ
+
η
u
2
+
δ
h
)
i
,
α
h
e
=
(
μ
h
+
ϕ
+
η
u
2
+
δ
h
)
i
,
e
=
(
μ
h
+
ϕ
+
η
u
2
+
δ
h
)
i
α
h
.
Equating these two equations for e we obtain
(27)
(
1
-
u
1
)
(
1
-
i
-
r
)
β
v
h
ϑ
z
+
κ
1
θ
μ
h
+
α
h
+
(
1
-
u
1
)
β
v
h
ϑ
z
=
(
μ
h
+
ϕ
+
η
u
2
+
δ
h
)
i
α
h
,
(
1
-
u
1
)
β
v
h
ϑ
α
h
z
i
+
[
μ
h
+
α
h
+
(
1
-
u
1
)
β
v
h
ϑ
z
(
μ
h
+
ϕ
+
η
u
2
+
δ
h
)
]
i
=
(
1
-
u
1
)
(
1
-
r
)
β
v
h
ϑ
α
h
z
+
κ
1
θ
α
h
+
μ
h
+
α
h
+
(
1
-
u
1
)
β
v
h
ϑ
z
.
Therefore, i with r and z becomes
(28)
i
=
(
(
1
-
u
1
)
(
1
-
r
)
β
v
h
ϑ
α
h
z
+
κ
1
θ
α
h
+
μ
h
+
α
h
+
(
1
-
u
1
)
β
v
h
ϑ
z
)
×
(
(
1
-
u
1
)
β
v
h
ϑ
α
h
z
i
+
[
μ
h
+
α
h
+
(
1
-
u
1
)
β
v
h
ϑ
z
(
μ
h
+
ϕ
+
η
u
2
+
δ
h
)
]
(
1
-
u
1
)
β
v
h
ϑ
α
h
z
i
)
-
1
.
Now solve for r from the third equation
(29)
(
ϕ
+
η
u
2
)
ρ
i
=
(
μ
h
+
ψ
)
r
,
r
=
(
ϕ
+
η
u
2
)
i
μ
h
+
ψ
.
Let a
1
=
(
ϕ
+
η
u
2
)
/
(
μ
h
+
ψ
), b
1
=
1
-
u
1, c
=
β
v
h
ϑ
α
h, d
=
κ
1
θ
α
h, g
=
μ
h
+
α
h, h
=
β
v
h
ϑ, and k
=
μ
h
+
ϕ
+
η
u
2
+
δ
h. Then substituting r in the i equation we obtain
(30)
i
=
b
1
(
1
-
a
1
i
)
c
z
+
d
+
g
+
b
h
z
b
1
c
z
+
g
+
b
1
h
k
z
,
(
b
1
c
z
+
g
+
b
1
h
k
z
+
b
1
c
a
1
z
b
1
c
z
+
g
+
b
1
h
k
z
)
i
=
b
1
c
z
+
d
+
g
+
b
1
h
z
b
1
c
z
+
g
+
b
1
h
k
z
,
i
=
b
1
c
z
+
d
+
g
+
b
1
h
z
b
1
c
z
+
g
+
b
1
h
k
z
+
b
1
c
a
1
z
.
Using equation five of the system (24), we obtain
(31)
(
μ
v
+
α
v
+
τ
u
3
)
y
+
β
h
v
i
y
=
β
h
v
ϑ
(
1
-
z
)
i
,
y
=
(
1
-
z
)
β
h
v
ϑ
i
μ
v
+
α
v
+
τ
u
3
+
β
h
v
ϑ
i
.
Solving for z gives
(32)
α
v
y
=
(
μ
v
+
τ
u
3
)
z
,
z
=
α
v
y
μ
v
+
τ
u
3
.
Let m
2
=
α
v
β
h
v
ϑ, n
2
=
μ
v
+
τ
u
3, p
1
=
μ
v
+
α
v
+
τ
u
3, and q
1
=
β
h
v
ϑ. Substituting y in the z equation gives
(33)
z
=
(
1
-
z
)
α
v
β
h
v
ϑ
i
(
μ
v
+
τ
u
3
)
(
μ
v
+
α
v
+
τ
u
3
+
β
h
v
ϑ
i
)
=
m
2
i
n
2
(
p
1
+
q
1
i
)
-
m
2
i
z
n
2
(
p
1
+
q
1
i
)
=
m
2
i
n
2
p
1
+
n
2
q
1
i
+
m
2
i
.
Substituting z in the i equation gives
(34)
i
=
d
+
g
+
(
b
1
c
+
b
1
h
)
(
m
2
i
/
(
n
2
p
1
+
n
2
q
1
i
+
m
2
i
)
)
g
+
(
b
1
c
+
b
1
h
k
+
b
1
c
a
1
)
(
m
2
i
/
(
n
2
p
1
+
n
2
q
1
i
+
m
2
i
)
)
.
The existence of the endemic equilibrium in Φ, can be determined when i
≠
0. After some algebraic manipulations we obtain
(35)
A
i
2
+
B
i
+
C
=
0
,
where
(36)
A
=
g
n
2
q
1
+
g
m
2
+
b
1
c
m
2
+
b
1
h
k
m
2
+
b
1
c
a
1
m
2
=
(
μ
h
+
α
h
)
(
μ
v
+
τ
u
3
)
β
h
v
ϑ
+
(
μ
h
+
α
h
)
α
v
β
h
v
ϑ
+
(
1
-
u
1
)
β
v
h
β
h
v
α
h
α
v
ϑ
2
+
(
1
-
u
1
)
(
μ
h
+
ϕ
+
η
u
2
+
δ
h
)
β
v
h
β
h
v
α
v
ϑ
2
+
(
1
-
u
1
)
(
ϕ
+
η
u
2
μ
h
+
ψ
)
β
v
h
β
h
v
α
h
α
v
ϑ
2
>
0
,
C
=
-
d
n
2
p
1
-
n
2
p
1
=
-
(
μ
v
+
τ
u
3
)
(
μ
v
+
α
v
+
τ
u
3
)
[
(
1
-
u
1
)
]
<
0
,
B
=
d
n
2
q
1
+
d
m
2
+
g
n
2
q
1
+
g
m
2
+
b
1
h
m
2
+
g
n
2
p
1
-
b
1
c
m
2
=
(
μ
v
+
τ
u
3
)
κ
1
θ
α
h
β
h
v
ϑ
+
κ
1
θ
α
h
α
v
β
h
v
ϑ
+
(
μ
h
+
α
h
)
(
μ
v
+
τ
u
3
)
β
h
v
ϑ
+
(
μ
h
+
α
h
)
α
v
β
h
v
ϑ
+
(
1
-
u
1
)
β
v
h
β
h
v
α
v
ϑ
2
+
(
μ
h
+
α
h
)
(
μ
v
+
τ
u
3
)
(
μ
v
+
α
v
+
τ
u
3
)
-
β
h
v
β
v
h
ϑ
2
α
h
α
v
(
1
-
u
1
)
=
(
μ
v
+
τ
u
3
)
κ
1
θ
α
h
β
h
v
ϑ
+
κ
1
θ
α
h
α
v
β
h
v
ϑ
+
(
μ
h
+
α
h
)
(
μ
v
+
τ
u
3
)
β
h
v
ϑ
+
(
μ
h
+
α
h
)
α
v
β
h
v
ϑ
+
(
1
-
u
1
)
β
v
h
β
h
v
α
v
ϑ
2
+
(
μ
h
+
α
h
)
(
μ
v
+
τ
u
3
)
(
μ
v
+
α
v
+
τ
u
3
)
(
1
-
R
e
2
)
.
For R
e
>
1, the existence of endemic equilibria is determined by the presence of positive real solutions of the quadratic expression (35). Consider
(37)
B
=
(
μ
v
+
τ
u
3
)
κ
1
θ
α
h
β
h
v
ϑ
+
κ
1
θ
α
h
α
v
β
h
v
ϑ
+
(
μ
h
+
α
h
)
(
μ
v
+
τ
u
3
)
β
h
v
ϑ
+
(
μ
h
+
α
h
)
κ
2
θ
α
v
β
h
v
ϑ
+
(
1
-
u
1
)
β
v
h
β
h
v
α
v
ϑ
2
+
(
μ
h
+
α
h
)
(
μ
v
+
τ
u
3
)
×
(
μ
v
+
α
v
+
τ
u
3
)
(
1
-
R
e
2
)
<
0
.
Since C
<
0 and A
>
0, then C
/
A
<
0. Thus, there exists exactly one positive endemic equilibrium for i
∈
(
0,1
] whenever R
e
>
1. This gives the threshold for the endemic persistence. Therefore, we have proved the existence and uniqueness of the endemic equilibrium, E
1, for the system (1). This result is summarized in the following theorem.Theorem 3.
IfR
e
>
1, then the model system (23) has a unique endemic equilibrium E
1.The result in Theorem3 indicates the impossibility of backward bifurcation in the malaria model, since it has no endemic equilibrium when R
e
<
1. Thus, the model (1) has a globally asymptotically stable disease-free equilibrium whenever R
e
≤
1.Next we need to investigate the global stability property of the endemic equilibrium for the case when the disease does not induce death.


#### 2.1.3. Global Stability of the Endemic Equilibrium forδ
h
=
0

The model system (1), when δ
h
=
0, has a unique endemic equilibrium. By letting
(38)
Φ
2
=
{
(
S
h
,
E
h
,
I
h
,
R
h
,
S
v
,
E
v
,
I
v
)
∈
R
+
7
:
E
h
=
I
h
l
l
l
l
=
E
v
=
I
v
=
0
(
S
h
,
E
h
,
I
h
,
R
h
,
S
v
,
E
v
,
I
v
)
∈
R
+
7
}
,
with
R
e
0
=
R
e
∣
δ
h
=
0
,
then the following can be claimed.Theorem 4.
The endemic equilibrium point of the malaria model (1) with δ
h
=
0 is globally asymptotically stable whenever R
e
0
>
1 in Φ and Φ
2.Proof.
As for the case of Theorem3, it can be shown that the unique endemic equilibrium for this special case exists only if R
e
0
>
1. Additionally, N
h
=
Λ
h
/
μ
h as t
→
∞. Letting S
h
=
Λ
h
/
μ
h
-
E
h
-
I
h
-
R
h and S
v
=
Λ
v
/
μ
v
-
E
v
-
I
v and substituting into (1) give the limiting system
(39)
d
E
h
d
t
=
(
1
-
u
1
)
λ
h
(
Λ
h
μ
h
-
E
h
-
I
h
-
R
h
)
+
κ
1
θ
Λ
h
μ
h
-
(
α
h
+
μ
h
)
E
h
,
d
I
h
d
t
=
α
h
E
h
-
(
ϕ
+
η
u
2
+
μ
h
+
δ
h
)
I
h
,
d
E
v
d
t
=
λ
v
(
Λ
v
μ
v
-
E
v
-
I
v
)
-
(
α
v
+
μ
v
+
τ
u
3
)
E
v
,
d
I
v
d
t
=
α
v
E
v
-
(
μ
v
+
τ
u
3
)
I
v
.
Dulac’s multiplier 1
/
E
h
I
h
E
v
I
v (see [19]) gives
(40)
∂
∂
E
h
[
(
1
-
u
1
)
β
v
h
ϑ
E
h
I
h
E
v
Λ
h
/
μ
h
(
Λ
h
μ
h
-
E
h
-
I
h
-
R
h
)
+
κ
1
θ
Λ
h
E
h
I
h
E
v
I
v
μ
h
-
(
α
h
+
μ
h
)
I
h
E
v
I
v
]
+
∂
∂
I
h
[
α
h
I
h
E
v
I
v
-
(
ϕ
+
η
u
2
+
μ
h
+
δ
h
)
E
h
E
v
I
v
]
+
∂
∂
E
v
[
(
α
v
+
μ
v
+
τ
u
3
)
E
h
I
h
E
v
β
h
v
ϑ
E
h
E
v
I
v
Λ
h
/
μ
h
(
Λ
v
μ
v
-
E
v
-
I
v
)
ddddddldd
-
(
α
v
+
μ
v
+
τ
u
3
)
E
h
I
h
E
v
]
+
∂
∂
I
v
[
α
v
E
h
I
h
I
v
-
(
μ
v
+
τ
u
3
)
E
h
I
h
E
v
]
=
-
β
v
h
ϑ
E
h
2
I
h
E
v
-
β
v
h
ϑ
μ
h
E
h
2
E
v
(
1
-
R
h
Λ
h
/
μ
h
)
-
κ
1
θ
Λ
h
E
h
2
I
h
E
v
I
v
μ
h
-
α
h
I
h
2
E
v
I
v
-
β
h
v
ϑ
E
h
E
v
2
I
v
[
1
-
μ
h
Λ
h
]
-
α
v
E
h
I
h
I
v
2
<
0
,
since E
h
+
I
h
+
R
h
<
Λ
h
/
μ
h and E
v
+
I
v
<
Λ
v
/
μ
v in Φ
2.
Thus, by Dulac’s criterion, there are no periodic orbits inΦ and Φ
2. Since Φ
2 is positively invariant and the endemic equilibrium exists whenever R
e
0
>
1, then from the Poincare-Bendixson Theorem [20], it follows that all solutions of the limiting system originating in Φ remain in Φ for all t. Furthermore, the absence of periodic orbits in Φ implies that the unique endemic equilibrium of the special case of the malaria model is globally asymptotically stable whenever R
e
0
>
1.The malaria model has a locally asymptotically stable disease-free equilibrium wheneverR
e
<
1 and a unique endemic equilibrium whenever R
e
>
1. In addition, the unique endemic equilibrium is globally asymptotically stable for the case δ
h
=
0 if R
e
0
>
1.In the next section, we apply the optimal control method using Pontryagin’s Maximum Principle to determine the necessary conditions for the combined optimal control of ITNs, IRS, and treatment effort which are being practised in Karonga District, Malawi.


## 2.1. Existence and Stability of Equilibrium Points

We analyse system (1) to obtain the equilibrium points of the system and their stability. Let (
S
h
*
,
E
h
*
,
I
h
*
,
R
h
*
,
S
v
*
,
E
v
*
,
I
v
*
) be the equilibrium points of system (1). At an equilibrium point, we have
(7)
S
h
′
(
t
)
=
E
h
′
(
t
)
=
I
h
′
(
t
)
=
R
h
′
(
t
)
=
S
v
′
(
t
)
=
E
v
′
(
t
)
=
I
v
′
(
t
)
=
0
.


### 2.1.1. Disease-Free Equilibrium,E
0

In the absence of malaria, that is,E
h
*
=
I
h
*
=
R
h
*
=
E
v
*
=
I
v
*
=
0, model system (1) has an equilibrium point called the disease-free equilibrium, E
0, and is given by
(8)
E
0
=
(
Λ
h
μ
h
,
0,0
,
0
,
Λ
h
μ
v
,
0,0
)
.
To establish the linear stability of E
0, we employ van den Driessche and Watmough’s next generation matrix approach [17]. A reproduction number obtained this way determines the local stability of the disease-free equilibrium point for R
e
<
1 and instability for R
e
>
1. Following van den Driessche and Watmough [17], the associated next generation matrices F and V of system (1) can be determined from F
i and V
i, respectively, where
(9)
F
i
=
[
(
1
-
u
1
)
β
v
h
ϑ
I
v
S
h
N
h
+
κ
1
θ
N
h
0
β
h
v
ϑ
I
h
S
v
N
h
0
]
,
V
i
=
[
(
α
h
+
μ
h
)
E
h
(
ϕ
+
η
u
2
+
μ
h
+
δ
h
)
I
h
-
α
h
E
h
(
α
v
+
μ
v
+
τ
u
3
)
E
v
(
μ
v
+
τ
u
3
)
I
v
-
α
v
E
v
]
.
The Jacobian matrix of F
i and V
i is given by
(10)
F
=
[
0
0
0
r
1
0
0
0
0
0
j
0
0
0
0
0
0
]
,
V
=
[
m
0
0
0
-
α
h
n
0
0
0
0
a
0
0
0
-
α
v
b
]
,
where
(11)
r
1
=
(
1
-
u
1
)
β
v
h
ϑ
,
j
=
β
h
v
ϑ
Λ
v
μ
h
Λ
h
μ
v
,
m
=
α
h
+
μ
h
,
n
=
ϕ
+
η
u
2
+
μ
h
+
δ
h
,
a
=
α
v
+
μ
v
+
τ
u
3
,
b
=
μ
v
+
τ
u
3
.
Algebraic manipulation of the matrices leads to the effective reproduction number
(12)
R
e
=
(
×
(
α
v
+
μ
v
+
τ
u
3
)
(
μ
v
+
τ
u
3
)
μ
v
Λ
h
ϑ
2
)
-
1
β
h
v
β
v
h
ϑ
2
α
h
α
v
(
1
-
u
1
)
Λ
v
μ
h
f
f
i
×
(
ϑ
2
(
α
h
+
μ
h
)
(
ϕ
+
η
u
2
+
μ
h
+
δ
h
)
f
f
f
i
f
i
×
(
α
v
+
μ
v
+
τ
u
3
)
(
μ
v
+
τ
u
3
)
μ
v
Λ
h
ϑ
2
)
-
1
)
1
/
2
,
where
(13)
R
e
v
=
β
h
v
α
v
ϑ
Λ
v
μ
v
(
α
v
+
μ
v
+
τ
u
3
)
(
μ
v
+
τ
u
3
)
is the contribution of the mosquito population when it infects the humans, and
(14)
R
e
h
=
β
v
h
ϑ
μ
h
α
h
(
1
-
u
1
)
Λ
h
(
α
h
+
μ
h
)
(
ϕ
+
η
u
2
+
μ
h
+
δ
h
)
is the human contribution when they infect the mosquitoes.The expression for the effective reproduction number,R
e, has a biological meaning that is readily interpreted from terms under the square root sign. Consider the following terms.(i)
β
v
h
ϑ
α
v
/
(
α
v
+
μ
v
+
τ
u
3
)
(
μ
v
+
τ
u
3
) represents the number of secondary human infections caused by one infected mosquito vector.
(ii)
β
h
v
ϑ
α
h
/
(
α
h
+
μ
h
)
(
ϕ
+
η
u
2
+
μ
h
+
δ
h
) represents the number of secondary mosquito infections caused by one infected human host.The square root represents the geometric mean of the average number of secondary host infections produced by one vector and the average number of secondary vector infections produced by one host. This effective reproduction number serves as an invasion threshold both for predicting outbreaks and evaluating control strategies that would reduce the spread of the disease. The threshold quantity,R
e, measures the average number of secondary cases generated by a single infected individual in a susceptible human population [18], where a fraction of the susceptible human population is under prevention and the infected class is under treatment. In the absence of any protective measure, the effective reproduction number R
e with treatment is
(15)
R
e
t
=
β
h
v
β
v
h
ϑ
2
α
h
α
v
Λ
v
μ
h
(
α
h
+
μ
h
)
(
ϕ
+
η
u
2
+
μ
h
+
δ
h
)
(
α
v
+
μ
v
)
μ
v
2
Λ
h
.
Also if ITNs are the only intervention strategy, then
(16)
R
e
N
=
β
h
v
β
v
h
ϑ
2
α
h
α
v
Λ
v
μ
h
(
1
-
u
1
)
(
α
h
+
μ
h
)
(
ϕ
+
μ
h
+
δ
h
)
(
α
v
+
μ
v
)
μ
v
2
Λ
h
.
Similarly, if IRS is the only means of protection, then
(17)
R
e
s
=
(
β
h
v
β
v
h
ϑ
2
α
h
α
v
μ
h
Λ
v
fffif
×
(
(
α
h
+
μ
h
)
(
ϕ
+
μ
h
+
δ
h
)
(
α
v
+
μ
v
+
τ
u
3
)
f
f
f
f
f
f
×
(
μ
v
+
τ
u
3
)
μ
v
Λ
h
)
-
1
)
1
/
2
.
The value of the basic reproduction number can be obtained from the value of the effective reproduction number when there are no control measures (
u
1
=
u
2
=
u
3
=
0
):
(18)
R
0
=
β
h
v
β
u
h
ϑ
2
α
h
α
v
Λ
v
μ
h
(
α
h
+
μ
h
)
(
ϕ
+
μ
h
+
δ
h
)
(
α
v
+
μ
v
)
μ
v
2
Λ
h
.
In general, it is easy to prove that
(19)
R
e
≤
R
0
for 0
≤
u
1, u
3
≤
1, due to reduction of likelihood of infection by protection. This implies that ITNs and IRS have a positive impact on the malaria dynamics as they contribute to the reduction of secondary infections. Therefore, from van den Driessche and Watmough [17] (Theorem 2), the following result holds.Lemma 2.
The disease-free equilibrium,E
0, of the malaria model with intervention strategies (1), given by (8), is locally asymptotically stable if R
e
<
1 and unstable if R
e
>
1.


### 2.1.2. Endemic Equilibrium Point,E
1

When malaria is present, the model system (1) has a steady state, E
1, called the endemic equilibrium. In order to establish the stability of E
1, we express system (1) in dimensionless variables, namely, S
h
/
N
h
=
s, E
h
/
N
h
=
e, I
h
/
N
h
=
i, R
h
/
N
h
=
r, S
v
/
N
v
=
x, E
v
/
N
v
=
y, and I
v
/
N
v
=
z, where d
N
h
/
d
t
=
Λ
h
-
μ
h
N
h
-
δ
h
I
h, d
N
v
/
d
t
=
Λ
v
-
μ
v
N
v. Then differentiating with respect to time, t, respectively results in the following system:
(20)
d
s
d
t
=
1
N
h
[
d
S
h
d
t
-
s
d
N
h
d
t
]
=
1
N
h
[
Λ
h
+
(
1
-
κ
1
)
θ
N
h
+
(
ϕ
+
η
u
2
)
(
1
-
ρ
)
i
N
h
f
f
f
f
i
f
-
(
1
-
u
1
)
β
v
h
ϑ
i
s
N
h
-
μ
h
s
N
h
+
ψ
r
N
h
]
.
Let m
1
=
(
1
-
κ
1
)
θ, n
1
=
(
ϕ
+
η
u
2
)
(
1
-
ρ
), and q
=
(
1
-
u
1
)
β
v
h
ϑ.Then(21)
d
s
d
t
=
1
N
h
[
Λ
h
+
m
1
N
h
+
n
1
i
N
h
-
q
i
s
N
h
-
u
h
s
N
h
+
ψ
r
N
h
]
-
s
N
h
[
Λ
h
-
δ
h
i
N
h
-
(
μ
h
-
θ
)
N
h
]
=
Λ
h
N
h
-
[
Λ
h
N
h
+
θ
-
δ
h
i
]
s
+
(
1
-
κ
1
)
θ
+
(
ϕ
+
η
u
2
)
(
1
-
ρ
)
i
-
(
1
-
u
1
)
β
v
h
ϑ
z
s
+
ψ
r
,
d
e
d
t
=
1
N
h
[
d
E
d
t
-
e
d
N
h
d
t
]
=
1
N
h
[
(
1
-
u
1
)
β
v
h
ϑ
z
s
N
h
+
κ
1
θ
N
h
-
α
h
e
N
h
-
μ
h
e
N
h
]
-
e
N
h
[
Λ
h
-
δ
h
i
N
h
-
(
μ
h
-
θ
)
N
h
]
=
(
1
-
u
1
)
β
v
h
ϑ
z
s
+
κ
1
θ
-
[
Λ
h
N
h
-
δ
h
i
+
α
h
+
θ
]
e
,
d
i
d
t
=
1
N
h
[
d
I
d
t
-
i
d
N
h
d
t
]
=
1
N
h
[
α
h
e
N
h
-
(
ϕ
+
η
u
2
)
(
1
-
ρ
)
i
N
h
fffffffffi
-
(
ϕ
+
η
u
2
)
ρ
i
N
h
-
(
μ
h
+
δ
h
)
i
N
h
]
-
i
N
h
[
Λ
h
-
δ
h
i
N
h
-
(
μ
h
-
θ
)
N
h
]
=
α
h
e
-
[
Λ
h
N
h
+
ϕ
+
η
u
2
+
δ
h
+
θ
-
δ
h
i
]
i
,
d
r
d
t
=
1
N
h
[
d
R
d
t
-
r
d
N
h
d
t
]
=
1
N
h
[
(
ϕ
+
η
u
2
)
ρ
i
N
h
-
(
μ
h
+
ψ
)
r
N
h
]
-
r
N
h
[
Λ
h
-
δ
h
i
N
h
-
(
μ
h
-
θ
)
N
h
]
=
(
ϕ
+
η
u
2
)
ρ
i
-
[
Λ
h
N
h
+
ψ
+
θ
-
δ
h
i
]
r
,
d
x
d
t
=
1
N
v
[
d
S
v
d
t
-
x
d
N
v
d
t
]
=
1
N
v
[
Λ
v
-
β
h
v
ϑ
i
x
N
v
-
(
μ
v
+
τ
u
3
)
x
N
v
]
-
x
N
v
[
Λ
v
-
μ
v
N
v
]
=
[
Λ
v
N
v
-
x
Λ
v
N
v
]
-
β
h
v
ϑ
i
x
-
τ
u
3
x
,
d
y
d
t
=
1
N
v
[
d
E
v
d
t
-
y
d
N
v
d
t
]
=
1
N
v
[
β
h
v
ϑ
i
x
N
v
-
α
v
y
N
v
-
(
μ
v
+
τ
u
3
)
y
N
v
]
-
y
N
v
[
Λ
v
-
μ
v
N
v
]
=
β
h
v
ϑ
i
x
-
[
Λ
v
N
v
+
α
v
+
τ
u
3
]
y
,
d
z
d
t
=
1
N
v
[
d
I
v
d
t
-
Z
d
N
v
d
t
]
=
1
N
v
[
α
v
y
N
v
-
(
μ
v
+
τ
u
3
)
z
N
v
]
-
z
N
v
[
Λ
v
-
μ
v
N
v
]
=
α
v
y
-
[
Λ
v
N
v
+
τ
u
3
]
z
,
d
N
h
d
t
=
[
Λ
h
N
h
-
μ
h
-
δ
h
i
]
N
h
,
d
N
v
d
t
=
[
Λ
v
N
v
-
μ
v
]
N
v
.
The system can now be reduced to a nine-dimensional system by eliminating s and x since s
=
1
-
e
-
i
-
r and x
=
1
-
y
-
z. The following feasible region
(22)
Φ
1
=
{
Λ
v
μ
v
e
,
i
,
r
,
N
h
,
y
,
z
,
N
v
∈
R
+
7
∣
e
≥
0
,
i
≥
0
,
r
≥
0
,
e
+
i
+
r
≤
1
,
N
h
≤
Λ
h
μ
h
,
y
≥
0
,
z
≥
0
,
y
+
z
≤
1
,
N
v
≤
Λ
v
μ
v
}
can be shown to be positively invariant. R
+
7 denotes the nonnegative cone of R
7 including its lower dimensional faces. Thus we have the following system of equations:
(23)
d
e
d
t
=
(
1
-
u
1
)
(
1
-
e
-
i
-
r
)
β
v
h
ϑ
z
+
κ
1
θ
-
[
Λ
h
N
h
-
δ
h
i
+
α
h
+
θ
]
e
,
d
i
d
t
=
α
h
e
-
[
Λ
h
N
h
+
ϕ
+
η
u
2
+
δ
h
+
θ
-
δ
h
i
]
i
,
d
r
d
t
=
(
ϕ
+
η
u
2
)
ρ
i
-
[
Λ
h
N
h
+
ψ
+
θ
+
-
δ
h
i
]
r
,
d
N
h
d
t
=
Λ
h
-
δ
h
i
N
h
-
μ
h
N
h
,
d
y
d
t
=
(
1
-
y
-
z
)
β
h
v
ϑ
i
-
[
Λ
v
N
v
+
α
v
+
τ
u
3
]
y
,
d
z
d
t
=
α
v
-
[
Λ
v
N
v
+
τ
u
3
]
z
,
d
N
v
d
t
=
Λ
v
-
μ
v
N
v
.
We seek to establish whether a unique endemic equilibrium exists. This is done by making more realistic assumptions, namely, that the protective control measures may not be totally effective.Existence and Uniqueness of Endemic Equilibrium,
E
1. To compute the steady states of the system (23), we set the derivative with respect to time in (23) equal to zero, and after simplification the following algebraic equations are obtained:
(24)
(
1
-
u
1
)
(
1
-
e
-
i
-
r
)
β
v
h
ϑ
z
+
κ
1
θ
=
[
Λ
h
N
h
-
δ
h
i
+
α
h
+
θ
]
e
,
α
h
e
=
[
Λ
h
N
h
+
ϕ
+
η
u
2
+
δ
h
+
θ
-
δ
h
i
]
i
,
(
ϕ
+
η
u
2
)
ρ
i
=
[
Λ
h
N
h
+
ψ
+
θ
-
δ
h
i
]
r
,
Λ
h
N
h
=
δ
h
i
+
μ
h
,
(
1
-
y
-
z
)
β
h
v
ϑ
i
=
[
Λ
v
N
v
+
α
v
+
τ
u
3
]
y
,
α
v
y
=
[
Λ
v
N
v
+
η
u
3
]
z
,
Λ
v
N
v
=
μ
v
.
To calculate the dimensionless proportions in terms of i, consider the first equation in the system (24):
(25)
(
1
-
i
-
r
)
p
z
-
p
z
e
+
κ
1
θ
=
[
μ
h
+
α
h
]
e
,
(
μ
h
+
α
h
)
e
+
(
1
-
u
1
)
β
v
h
ϑ
z
e
=
(
1
-
u
1
)
(
1
-
i
-
r
)
β
v
h
ϑ
z
+
κ
1
θ
,
(
μ
h
+
α
h
+
(
1
-
u
1
)
β
v
h
z
)
e
=
(
1
-
u
1
)
(
1
-
i
-
r
)
β
v
h
ϑ
z
+
κ
1
θ
,
e
=
(
1
-
u
1
)
(
1
-
i
-
r
)
β
v
h
ϑ
z
+
κ
1
θ
μ
h
+
α
h
+
(
1
-
u
1
)
β
v
h
ϑ
z
,
where p
=
(
1
-
u
1
)
β
v
h
ϑ.From the second equation we have(26)
α
h
e
=
(
μ
h
+
ϕ
+
η
u
2
+
δ
h
)
i
,
α
h
e
=
(
μ
h
+
ϕ
+
η
u
2
+
δ
h
)
i
,
e
=
(
μ
h
+
ϕ
+
η
u
2
+
δ
h
)
i
α
h
.
Equating these two equations for e we obtain
(27)
(
1
-
u
1
)
(
1
-
i
-
r
)
β
v
h
ϑ
z
+
κ
1
θ
μ
h
+
α
h
+
(
1
-
u
1
)
β
v
h
ϑ
z
=
(
μ
h
+
ϕ
+
η
u
2
+
δ
h
)
i
α
h
,
(
1
-
u
1
)
β
v
h
ϑ
α
h
z
i
+
[
μ
h
+
α
h
+
(
1
-
u
1
)
β
v
h
ϑ
z
(
μ
h
+
ϕ
+
η
u
2
+
δ
h
)
]
i
=
(
1
-
u
1
)
(
1
-
r
)
β
v
h
ϑ
α
h
z
+
κ
1
θ
α
h
+
μ
h
+
α
h
+
(
1
-
u
1
)
β
v
h
ϑ
z
.
Therefore, i with r and z becomes
(28)
i
=
(
(
1
-
u
1
)
(
1
-
r
)
β
v
h
ϑ
α
h
z
+
κ
1
θ
α
h
+
μ
h
+
α
h
+
(
1
-
u
1
)
β
v
h
ϑ
z
)
×
(
(
1
-
u
1
)
β
v
h
ϑ
α
h
z
i
+
[
μ
h
+
α
h
+
(
1
-
u
1
)
β
v
h
ϑ
z
(
μ
h
+
ϕ
+
η
u
2
+
δ
h
)
]
(
1
-
u
1
)
β
v
h
ϑ
α
h
z
i
)
-
1
.
Now solve for r from the third equation
(29)
(
ϕ
+
η
u
2
)
ρ
i
=
(
μ
h
+
ψ
)
r
,
r
=
(
ϕ
+
η
u
2
)
i
μ
h
+
ψ
.
Let a
1
=
(
ϕ
+
η
u
2
)
/
(
μ
h
+
ψ
), b
1
=
1
-
u
1, c
=
β
v
h
ϑ
α
h, d
=
κ
1
θ
α
h, g
=
μ
h
+
α
h, h
=
β
v
h
ϑ, and k
=
μ
h
+
ϕ
+
η
u
2
+
δ
h. Then substituting r in the i equation we obtain
(30)
i
=
b
1
(
1
-
a
1
i
)
c
z
+
d
+
g
+
b
h
z
b
1
c
z
+
g
+
b
1
h
k
z
,
(
b
1
c
z
+
g
+
b
1
h
k
z
+
b
1
c
a
1
z
b
1
c
z
+
g
+
b
1
h
k
z
)
i
=
b
1
c
z
+
d
+
g
+
b
1
h
z
b
1
c
z
+
g
+
b
1
h
k
z
,
i
=
b
1
c
z
+
d
+
g
+
b
1
h
z
b
1
c
z
+
g
+
b
1
h
k
z
+
b
1
c
a
1
z
.
Using equation five of the system (24), we obtain
(31)
(
μ
v
+
α
v
+
τ
u
3
)
y
+
β
h
v
i
y
=
β
h
v
ϑ
(
1
-
z
)
i
,
y
=
(
1
-
z
)
β
h
v
ϑ
i
μ
v
+
α
v
+
τ
u
3
+
β
h
v
ϑ
i
.
Solving for z gives
(32)
α
v
y
=
(
μ
v
+
τ
u
3
)
z
,
z
=
α
v
y
μ
v
+
τ
u
3
.
Let m
2
=
α
v
β
h
v
ϑ, n
2
=
μ
v
+
τ
u
3, p
1
=
μ
v
+
α
v
+
τ
u
3, and q
1
=
β
h
v
ϑ. Substituting y in the z equation gives
(33)
z
=
(
1
-
z
)
α
v
β
h
v
ϑ
i
(
μ
v
+
τ
u
3
)
(
μ
v
+
α
v
+
τ
u
3
+
β
h
v
ϑ
i
)
=
m
2
i
n
2
(
p
1
+
q
1
i
)
-
m
2
i
z
n
2
(
p
1
+
q
1
i
)
=
m
2
i
n
2
p
1
+
n
2
q
1
i
+
m
2
i
.
Substituting z in the i equation gives
(34)
i
=
d
+
g
+
(
b
1
c
+
b
1
h
)
(
m
2
i
/
(
n
2
p
1
+
n
2
q
1
i
+
m
2
i
)
)
g
+
(
b
1
c
+
b
1
h
k
+
b
1
c
a
1
)
(
m
2
i
/
(
n
2
p
1
+
n
2
q
1
i
+
m
2
i
)
)
.
The existence of the endemic equilibrium in Φ, can be determined when i
≠
0. After some algebraic manipulations we obtain
(35)
A
i
2
+
B
i
+
C
=
0
,
where
(36)
A
=
g
n
2
q
1
+
g
m
2
+
b
1
c
m
2
+
b
1
h
k
m
2
+
b
1
c
a
1
m
2
=
(
μ
h
+
α
h
)
(
μ
v
+
τ
u
3
)
β
h
v
ϑ
+
(
μ
h
+
α
h
)
α
v
β
h
v
ϑ
+
(
1
-
u
1
)
β
v
h
β
h
v
α
h
α
v
ϑ
2
+
(
1
-
u
1
)
(
μ
h
+
ϕ
+
η
u
2
+
δ
h
)
β
v
h
β
h
v
α
v
ϑ
2
+
(
1
-
u
1
)
(
ϕ
+
η
u
2
μ
h
+
ψ
)
β
v
h
β
h
v
α
h
α
v
ϑ
2
>
0
,
C
=
-
d
n
2
p
1
-
n
2
p
1
=
-
(
μ
v
+
τ
u
3
)
(
μ
v
+
α
v
+
τ
u
3
)
[
(
1
-
u
1
)
]
<
0
,
B
=
d
n
2
q
1
+
d
m
2
+
g
n
2
q
1
+
g
m
2
+
b
1
h
m
2
+
g
n
2
p
1
-
b
1
c
m
2
=
(
μ
v
+
τ
u
3
)
κ
1
θ
α
h
β
h
v
ϑ
+
κ
1
θ
α
h
α
v
β
h
v
ϑ
+
(
μ
h
+
α
h
)
(
μ
v
+
τ
u
3
)
β
h
v
ϑ
+
(
μ
h
+
α
h
)
α
v
β
h
v
ϑ
+
(
1
-
u
1
)
β
v
h
β
h
v
α
v
ϑ
2
+
(
μ
h
+
α
h
)
(
μ
v
+
τ
u
3
)
(
μ
v
+
α
v
+
τ
u
3
)
-
β
h
v
β
v
h
ϑ
2
α
h
α
v
(
1
-
u
1
)
=
(
μ
v
+
τ
u
3
)
κ
1
θ
α
h
β
h
v
ϑ
+
κ
1
θ
α
h
α
v
β
h
v
ϑ
+
(
μ
h
+
α
h
)
(
μ
v
+
τ
u
3
)
β
h
v
ϑ
+
(
μ
h
+
α
h
)
α
v
β
h
v
ϑ
+
(
1
-
u
1
)
β
v
h
β
h
v
α
v
ϑ
2
+
(
μ
h
+
α
h
)
(
μ
v
+
τ
u
3
)
(
μ
v
+
α
v
+
τ
u
3
)
(
1
-
R
e
2
)
.
For R
e
>
1, the existence of endemic equilibria is determined by the presence of positive real solutions of the quadratic expression (35). Consider
(37)
B
=
(
μ
v
+
τ
u
3
)
κ
1
θ
α
h
β
h
v
ϑ
+
κ
1
θ
α
h
α
v
β
h
v
ϑ
+
(
μ
h
+
α
h
)
(
μ
v
+
τ
u
3
)
β
h
v
ϑ
+
(
μ
h
+
α
h
)
κ
2
θ
α
v
β
h
v
ϑ
+
(
1
-
u
1
)
β
v
h
β
h
v
α
v
ϑ
2
+
(
μ
h
+
α
h
)
(
μ
v
+
τ
u
3
)
×
(
μ
v
+
α
v
+
τ
u
3
)
(
1
-
R
e
2
)
<
0
.
Since C
<
0 and A
>
0, then C
/
A
<
0. Thus, there exists exactly one positive endemic equilibrium for i
∈
(
0,1
] whenever R
e
>
1. This gives the threshold for the endemic persistence. Therefore, we have proved the existence and uniqueness of the endemic equilibrium, E
1, for the system (1). This result is summarized in the following theorem.Theorem 3.
IfR
e
>
1, then the model system (23) has a unique endemic equilibrium E
1.The result in Theorem3 indicates the impossibility of backward bifurcation in the malaria model, since it has no endemic equilibrium when R
e
<
1. Thus, the model (1) has a globally asymptotically stable disease-free equilibrium whenever R
e
≤
1.Next we need to investigate the global stability property of the endemic equilibrium for the case when the disease does not induce death.


### 2.1.3. Global Stability of the Endemic Equilibrium forδ
h
=
0

The model system (1), when δ
h
=
0, has a unique endemic equilibrium. By letting
(38)
Φ
2
=
{
(
S
h
,
E
h
,
I
h
,
R
h
,
S
v
,
E
v
,
I
v
)
∈
R
+
7
:
E
h
=
I
h
l
l
l
l
=
E
v
=
I
v
=
0
(
S
h
,
E
h
,
I
h
,
R
h
,
S
v
,
E
v
,
I
v
)
∈
R
+
7
}
,
with
R
e
0
=
R
e
∣
δ
h
=
0
,
then the following can be claimed.Theorem 4.
The endemic equilibrium point of the malaria model (1) with δ
h
=
0 is globally asymptotically stable whenever R
e
0
>
1 in Φ and Φ
2.Proof.
As for the case of Theorem3, it can be shown that the unique endemic equilibrium for this special case exists only if R
e
0
>
1. Additionally, N
h
=
Λ
h
/
μ
h as t
→
∞. Letting S
h
=
Λ
h
/
μ
h
-
E
h
-
I
h
-
R
h and S
v
=
Λ
v
/
μ
v
-
E
v
-
I
v and substituting into (1) give the limiting system
(39)
d
E
h
d
t
=
(
1
-
u
1
)
λ
h
(
Λ
h
μ
h
-
E
h
-
I
h
-
R
h
)
+
κ
1
θ
Λ
h
μ
h
-
(
α
h
+
μ
h
)
E
h
,
d
I
h
d
t
=
α
h
E
h
-
(
ϕ
+
η
u
2
+
μ
h
+
δ
h
)
I
h
,
d
E
v
d
t
=
λ
v
(
Λ
v
μ
v
-
E
v
-
I
v
)
-
(
α
v
+
μ
v
+
τ
u
3
)
E
v
,
d
I
v
d
t
=
α
v
E
v
-
(
μ
v
+
τ
u
3
)
I
v
.
Dulac’s multiplier 1
/
E
h
I
h
E
v
I
v (see [19]) gives
(40)
∂
∂
E
h
[
(
1
-
u
1
)
β
v
h
ϑ
E
h
I
h
E
v
Λ
h
/
μ
h
(
Λ
h
μ
h
-
E
h
-
I
h
-
R
h
)
+
κ
1
θ
Λ
h
E
h
I
h
E
v
I
v
μ
h
-
(
α
h
+
μ
h
)
I
h
E
v
I
v
]
+
∂
∂
I
h
[
α
h
I
h
E
v
I
v
-
(
ϕ
+
η
u
2
+
μ
h
+
δ
h
)
E
h
E
v
I
v
]
+
∂
∂
E
v
[
(
α
v
+
μ
v
+
τ
u
3
)
E
h
I
h
E
v
β
h
v
ϑ
E
h
E
v
I
v
Λ
h
/
μ
h
(
Λ
v
μ
v
-
E
v
-
I
v
)
ddddddldd
-
(
α
v
+
μ
v
+
τ
u
3
)
E
h
I
h
E
v
]
+
∂
∂
I
v
[
α
v
E
h
I
h
I
v
-
(
μ
v
+
τ
u
3
)
E
h
I
h
E
v
]
=
-
β
v
h
ϑ
E
h
2
I
h
E
v
-
β
v
h
ϑ
μ
h
E
h
2
E
v
(
1
-
R
h
Λ
h
/
μ
h
)
-
κ
1
θ
Λ
h
E
h
2
I
h
E
v
I
v
μ
h
-
α
h
I
h
2
E
v
I
v
-
β
h
v
ϑ
E
h
E
v
2
I
v
[
1
-
μ
h
Λ
h
]
-
α
v
E
h
I
h
I
v
2
<
0
,
since E
h
+
I
h
+
R
h
<
Λ
h
/
μ
h and E
v
+
I
v
<
Λ
v
/
μ
v in Φ
2.
Thus, by Dulac’s criterion, there are no periodic orbits inΦ and Φ
2. Since Φ
2 is positively invariant and the endemic equilibrium exists whenever R
e
0
>
1, then from the Poincare-Bendixson Theorem [20], it follows that all solutions of the limiting system originating in Φ remain in Φ for all t. Furthermore, the absence of periodic orbits in Φ implies that the unique endemic equilibrium of the special case of the malaria model is globally asymptotically stable whenever R
e
0
>
1.The malaria model has a locally asymptotically stable disease-free equilibrium wheneverR
e
<
1 and a unique endemic equilibrium whenever R
e
>
1. In addition, the unique endemic equilibrium is globally asymptotically stable for the case δ
h
=
0 if R
e
0
>
1.In the next section, we apply the optimal control method using Pontryagin’s Maximum Principle to determine the necessary conditions for the combined optimal control of ITNs, IRS, and treatment effort which are being practised in Karonga District, Malawi.


## 2.1.1. Disease-Free Equilibrium,E
0

In the absence of malaria, that is,E
h
*
=
I
h
*
=
R
h
*
=
E
v
*
=
I
v
*
=
0, model system (1) has an equilibrium point called the disease-free equilibrium, E
0, and is given by
(8)
E
0
=
(
Λ
h
μ
h
,
0,0
,
0
,
Λ
h
μ
v
,
0,0
)
.
To establish the linear stability of E
0, we employ van den Driessche and Watmough’s next generation matrix approach [17]. A reproduction number obtained this way determines the local stability of the disease-free equilibrium point for R
e
<
1 and instability for R
e
>
1. Following van den Driessche and Watmough [17], the associated next generation matrices F and V of system (1) can be determined from F
i and V
i, respectively, where
(9)
F
i
=
[
(
1
-
u
1
)
β
v
h
ϑ
I
v
S
h
N
h
+
κ
1
θ
N
h
0
β
h
v
ϑ
I
h
S
v
N
h
0
]
,
V
i
=
[
(
α
h
+
μ
h
)
E
h
(
ϕ
+
η
u
2
+
μ
h
+
δ
h
)
I
h
-
α
h
E
h
(
α
v
+
μ
v
+
τ
u
3
)
E
v
(
μ
v
+
τ
u
3
)
I
v
-
α
v
E
v
]
.
The Jacobian matrix of F
i and V
i is given by
(10)
F
=
[
0
0
0
r
1
0
0
0
0
0
j
0
0
0
0
0
0
]
,
V
=
[
m
0
0
0
-
α
h
n
0
0
0
0
a
0
0
0
-
α
v
b
]
,
where
(11)
r
1
=
(
1
-
u
1
)
β
v
h
ϑ
,
j
=
β
h
v
ϑ
Λ
v
μ
h
Λ
h
μ
v
,
m
=
α
h
+
μ
h
,
n
=
ϕ
+
η
u
2
+
μ
h
+
δ
h
,
a
=
α
v
+
μ
v
+
τ
u
3
,
b
=
μ
v
+
τ
u
3
.
Algebraic manipulation of the matrices leads to the effective reproduction number
(12)
R
e
=
(
×
(
α
v
+
μ
v
+
τ
u
3
)
(
μ
v
+
τ
u
3
)
μ
v
Λ
h
ϑ
2
)
-
1
β
h
v
β
v
h
ϑ
2
α
h
α
v
(
1
-
u
1
)
Λ
v
μ
h
f
f
i
×
(
ϑ
2
(
α
h
+
μ
h
)
(
ϕ
+
η
u
2
+
μ
h
+
δ
h
)
f
f
f
i
f
i
×
(
α
v
+
μ
v
+
τ
u
3
)
(
μ
v
+
τ
u
3
)
μ
v
Λ
h
ϑ
2
)
-
1
)
1
/
2
,
where
(13)
R
e
v
=
β
h
v
α
v
ϑ
Λ
v
μ
v
(
α
v
+
μ
v
+
τ
u
3
)
(
μ
v
+
τ
u
3
)
is the contribution of the mosquito population when it infects the humans, and
(14)
R
e
h
=
β
v
h
ϑ
μ
h
α
h
(
1
-
u
1
)
Λ
h
(
α
h
+
μ
h
)
(
ϕ
+
η
u
2
+
μ
h
+
δ
h
)
is the human contribution when they infect the mosquitoes.The expression for the effective reproduction number,R
e, has a biological meaning that is readily interpreted from terms under the square root sign. Consider the following terms.(i)
β
v
h
ϑ
α
v
/
(
α
v
+
μ
v
+
τ
u
3
)
(
μ
v
+
τ
u
3
) represents the number of secondary human infections caused by one infected mosquito vector.
(ii)
β
h
v
ϑ
α
h
/
(
α
h
+
μ
h
)
(
ϕ
+
η
u
2
+
μ
h
+
δ
h
) represents the number of secondary mosquito infections caused by one infected human host.The square root represents the geometric mean of the average number of secondary host infections produced by one vector and the average number of secondary vector infections produced by one host. This effective reproduction number serves as an invasion threshold both for predicting outbreaks and evaluating control strategies that would reduce the spread of the disease. The threshold quantity,R
e, measures the average number of secondary cases generated by a single infected individual in a susceptible human population [18], where a fraction of the susceptible human population is under prevention and the infected class is under treatment. In the absence of any protective measure, the effective reproduction number R
e with treatment is
(15)
R
e
t
=
β
h
v
β
v
h
ϑ
2
α
h
α
v
Λ
v
μ
h
(
α
h
+
μ
h
)
(
ϕ
+
η
u
2
+
μ
h
+
δ
h
)
(
α
v
+
μ
v
)
μ
v
2
Λ
h
.
Also if ITNs are the only intervention strategy, then
(16)
R
e
N
=
β
h
v
β
v
h
ϑ
2
α
h
α
v
Λ
v
μ
h
(
1
-
u
1
)
(
α
h
+
μ
h
)
(
ϕ
+
μ
h
+
δ
h
)
(
α
v
+
μ
v
)
μ
v
2
Λ
h
.
Similarly, if IRS is the only means of protection, then
(17)
R
e
s
=
(
β
h
v
β
v
h
ϑ
2
α
h
α
v
μ
h
Λ
v
fffif
×
(
(
α
h
+
μ
h
)
(
ϕ
+
μ
h
+
δ
h
)
(
α
v
+
μ
v
+
τ
u
3
)
f
f
f
f
f
f
×
(
μ
v
+
τ
u
3
)
μ
v
Λ
h
)
-
1
)
1
/
2
.
The value of the basic reproduction number can be obtained from the value of the effective reproduction number when there are no control measures (
u
1
=
u
2
=
u
3
=
0
):
(18)
R
0
=
β
h
v
β
u
h
ϑ
2
α
h
α
v
Λ
v
μ
h
(
α
h
+
μ
h
)
(
ϕ
+
μ
h
+
δ
h
)
(
α
v
+
μ
v
)
μ
v
2
Λ
h
.
In general, it is easy to prove that
(19)
R
e
≤
R
0
for 0
≤
u
1, u
3
≤
1, due to reduction of likelihood of infection by protection. This implies that ITNs and IRS have a positive impact on the malaria dynamics as they contribute to the reduction of secondary infections. Therefore, from van den Driessche and Watmough [17] (Theorem 2), the following result holds.Lemma 2.
The disease-free equilibrium,E
0, of the malaria model with intervention strategies (1), given by (8), is locally asymptotically stable if R
e
<
1 and unstable if R
e
>
1.


## 2.1.2. Endemic Equilibrium Point,E
1

When malaria is present, the model system (1) has a steady state, E
1, called the endemic equilibrium. In order to establish the stability of E
1, we express system (1) in dimensionless variables, namely, S
h
/
N
h
=
s, E
h
/
N
h
=
e, I
h
/
N
h
=
i, R
h
/
N
h
=
r, S
v
/
N
v
=
x, E
v
/
N
v
=
y, and I
v
/
N
v
=
z, where d
N
h
/
d
t
=
Λ
h
-
μ
h
N
h
-
δ
h
I
h, d
N
v
/
d
t
=
Λ
v
-
μ
v
N
v. Then differentiating with respect to time, t, respectively results in the following system:
(20)
d
s
d
t
=
1
N
h
[
d
S
h
d
t
-
s
d
N
h
d
t
]
=
1
N
h
[
Λ
h
+
(
1
-
κ
1
)
θ
N
h
+
(
ϕ
+
η
u
2
)
(
1
-
ρ
)
i
N
h
f
f
f
f
i
f
-
(
1
-
u
1
)
β
v
h
ϑ
i
s
N
h
-
μ
h
s
N
h
+
ψ
r
N
h
]
.
Let m
1
=
(
1
-
κ
1
)
θ, n
1
=
(
ϕ
+
η
u
2
)
(
1
-
ρ
), and q
=
(
1
-
u
1
)
β
v
h
ϑ.Then(21)
d
s
d
t
=
1
N
h
[
Λ
h
+
m
1
N
h
+
n
1
i
N
h
-
q
i
s
N
h
-
u
h
s
N
h
+
ψ
r
N
h
]
-
s
N
h
[
Λ
h
-
δ
h
i
N
h
-
(
μ
h
-
θ
)
N
h
]
=
Λ
h
N
h
-
[
Λ
h
N
h
+
θ
-
δ
h
i
]
s
+
(
1
-
κ
1
)
θ
+
(
ϕ
+
η
u
2
)
(
1
-
ρ
)
i
-
(
1
-
u
1
)
β
v
h
ϑ
z
s
+
ψ
r
,
d
e
d
t
=
1
N
h
[
d
E
d
t
-
e
d
N
h
d
t
]
=
1
N
h
[
(
1
-
u
1
)
β
v
h
ϑ
z
s
N
h
+
κ
1
θ
N
h
-
α
h
e
N
h
-
μ
h
e
N
h
]
-
e
N
h
[
Λ
h
-
δ
h
i
N
h
-
(
μ
h
-
θ
)
N
h
]
=
(
1
-
u
1
)
β
v
h
ϑ
z
s
+
κ
1
θ
-
[
Λ
h
N
h
-
δ
h
i
+
α
h
+
θ
]
e
,
d
i
d
t
=
1
N
h
[
d
I
d
t
-
i
d
N
h
d
t
]
=
1
N
h
[
α
h
e
N
h
-
(
ϕ
+
η
u
2
)
(
1
-
ρ
)
i
N
h
fffffffffi
-
(
ϕ
+
η
u
2
)
ρ
i
N
h
-
(
μ
h
+
δ
h
)
i
N
h
]
-
i
N
h
[
Λ
h
-
δ
h
i
N
h
-
(
μ
h
-
θ
)
N
h
]
=
α
h
e
-
[
Λ
h
N
h
+
ϕ
+
η
u
2
+
δ
h
+
θ
-
δ
h
i
]
i
,
d
r
d
t
=
1
N
h
[
d
R
d
t
-
r
d
N
h
d
t
]
=
1
N
h
[
(
ϕ
+
η
u
2
)
ρ
i
N
h
-
(
μ
h
+
ψ
)
r
N
h
]
-
r
N
h
[
Λ
h
-
δ
h
i
N
h
-
(
μ
h
-
θ
)
N
h
]
=
(
ϕ
+
η
u
2
)
ρ
i
-
[
Λ
h
N
h
+
ψ
+
θ
-
δ
h
i
]
r
,
d
x
d
t
=
1
N
v
[
d
S
v
d
t
-
x
d
N
v
d
t
]
=
1
N
v
[
Λ
v
-
β
h
v
ϑ
i
x
N
v
-
(
μ
v
+
τ
u
3
)
x
N
v
]
-
x
N
v
[
Λ
v
-
μ
v
N
v
]
=
[
Λ
v
N
v
-
x
Λ
v
N
v
]
-
β
h
v
ϑ
i
x
-
τ
u
3
x
,
d
y
d
t
=
1
N
v
[
d
E
v
d
t
-
y
d
N
v
d
t
]
=
1
N
v
[
β
h
v
ϑ
i
x
N
v
-
α
v
y
N
v
-
(
μ
v
+
τ
u
3
)
y
N
v
]
-
y
N
v
[
Λ
v
-
μ
v
N
v
]
=
β
h
v
ϑ
i
x
-
[
Λ
v
N
v
+
α
v
+
τ
u
3
]
y
,
d
z
d
t
=
1
N
v
[
d
I
v
d
t
-
Z
d
N
v
d
t
]
=
1
N
v
[
α
v
y
N
v
-
(
μ
v
+
τ
u
3
)
z
N
v
]
-
z
N
v
[
Λ
v
-
μ
v
N
v
]
=
α
v
y
-
[
Λ
v
N
v
+
τ
u
3
]
z
,
d
N
h
d
t
=
[
Λ
h
N
h
-
μ
h
-
δ
h
i
]
N
h
,
d
N
v
d
t
=
[
Λ
v
N
v
-
μ
v
]
N
v
.
The system can now be reduced to a nine-dimensional system by eliminating s and x since s
=
1
-
e
-
i
-
r and x
=
1
-
y
-
z. The following feasible region
(22)
Φ
1
=
{
Λ
v
μ
v
e
,
i
,
r
,
N
h
,
y
,
z
,
N
v
∈
R
+
7
∣
e
≥
0
,
i
≥
0
,
r
≥
0
,
e
+
i
+
r
≤
1
,
N
h
≤
Λ
h
μ
h
,
y
≥
0
,
z
≥
0
,
y
+
z
≤
1
,
N
v
≤
Λ
v
μ
v
}
can be shown to be positively invariant. R
+
7 denotes the nonnegative cone of R
7 including its lower dimensional faces. Thus we have the following system of equations:
(23)
d
e
d
t
=
(
1
-
u
1
)
(
1
-
e
-
i
-
r
)
β
v
h
ϑ
z
+
κ
1
θ
-
[
Λ
h
N
h
-
δ
h
i
+
α
h
+
θ
]
e
,
d
i
d
t
=
α
h
e
-
[
Λ
h
N
h
+
ϕ
+
η
u
2
+
δ
h
+
θ
-
δ
h
i
]
i
,
d
r
d
t
=
(
ϕ
+
η
u
2
)
ρ
i
-
[
Λ
h
N
h
+
ψ
+
θ
+
-
δ
h
i
]
r
,
d
N
h
d
t
=
Λ
h
-
δ
h
i
N
h
-
μ
h
N
h
,
d
y
d
t
=
(
1
-
y
-
z
)
β
h
v
ϑ
i
-
[
Λ
v
N
v
+
α
v
+
τ
u
3
]
y
,
d
z
d
t
=
α
v
-
[
Λ
v
N
v
+
τ
u
3
]
z
,
d
N
v
d
t
=
Λ
v
-
μ
v
N
v
.
We seek to establish whether a unique endemic equilibrium exists. This is done by making more realistic assumptions, namely, that the protective control measures may not be totally effective.Existence and Uniqueness of Endemic Equilibrium,
E
1. To compute the steady states of the system (23), we set the derivative with respect to time in (23) equal to zero, and after simplification the following algebraic equations are obtained:
(24)
(
1
-
u
1
)
(
1
-
e
-
i
-
r
)
β
v
h
ϑ
z
+
κ
1
θ
=
[
Λ
h
N
h
-
δ
h
i
+
α
h
+
θ
]
e
,
α
h
e
=
[
Λ
h
N
h
+
ϕ
+
η
u
2
+
δ
h
+
θ
-
δ
h
i
]
i
,
(
ϕ
+
η
u
2
)
ρ
i
=
[
Λ
h
N
h
+
ψ
+
θ
-
δ
h
i
]
r
,
Λ
h
N
h
=
δ
h
i
+
μ
h
,
(
1
-
y
-
z
)
β
h
v
ϑ
i
=
[
Λ
v
N
v
+
α
v
+
τ
u
3
]
y
,
α
v
y
=
[
Λ
v
N
v
+
η
u
3
]
z
,
Λ
v
N
v
=
μ
v
.
To calculate the dimensionless proportions in terms of i, consider the first equation in the system (24):
(25)
(
1
-
i
-
r
)
p
z
-
p
z
e
+
κ
1
θ
=
[
μ
h
+
α
h
]
e
,
(
μ
h
+
α
h
)
e
+
(
1
-
u
1
)
β
v
h
ϑ
z
e
=
(
1
-
u
1
)
(
1
-
i
-
r
)
β
v
h
ϑ
z
+
κ
1
θ
,
(
μ
h
+
α
h
+
(
1
-
u
1
)
β
v
h
z
)
e
=
(
1
-
u
1
)
(
1
-
i
-
r
)
β
v
h
ϑ
z
+
κ
1
θ
,
e
=
(
1
-
u
1
)
(
1
-
i
-
r
)
β
v
h
ϑ
z
+
κ
1
θ
μ
h
+
α
h
+
(
1
-
u
1
)
β
v
h
ϑ
z
,
where p
=
(
1
-
u
1
)
β
v
h
ϑ.From the second equation we have(26)
α
h
e
=
(
μ
h
+
ϕ
+
η
u
2
+
δ
h
)
i
,
α
h
e
=
(
μ
h
+
ϕ
+
η
u
2
+
δ
h
)
i
,
e
=
(
μ
h
+
ϕ
+
η
u
2
+
δ
h
)
i
α
h
.
Equating these two equations for e we obtain
(27)
(
1
-
u
1
)
(
1
-
i
-
r
)
β
v
h
ϑ
z
+
κ
1
θ
μ
h
+
α
h
+
(
1
-
u
1
)
β
v
h
ϑ
z
=
(
μ
h
+
ϕ
+
η
u
2
+
δ
h
)
i
α
h
,
(
1
-
u
1
)
β
v
h
ϑ
α
h
z
i
+
[
μ
h
+
α
h
+
(
1
-
u
1
)
β
v
h
ϑ
z
(
μ
h
+
ϕ
+
η
u
2
+
δ
h
)
]
i
=
(
1
-
u
1
)
(
1
-
r
)
β
v
h
ϑ
α
h
z
+
κ
1
θ
α
h
+
μ
h
+
α
h
+
(
1
-
u
1
)
β
v
h
ϑ
z
.
Therefore, i with r and z becomes
(28)
i
=
(
(
1
-
u
1
)
(
1
-
r
)
β
v
h
ϑ
α
h
z
+
κ
1
θ
α
h
+
μ
h
+
α
h
+
(
1
-
u
1
)
β
v
h
ϑ
z
)
×
(
(
1
-
u
1
)
β
v
h
ϑ
α
h
z
i
+
[
μ
h
+
α
h
+
(
1
-
u
1
)
β
v
h
ϑ
z
(
μ
h
+
ϕ
+
η
u
2
+
δ
h
)
]
(
1
-
u
1
)
β
v
h
ϑ
α
h
z
i
)
-
1
.
Now solve for r from the third equation
(29)
(
ϕ
+
η
u
2
)
ρ
i
=
(
μ
h
+
ψ
)
r
,
r
=
(
ϕ
+
η
u
2
)
i
μ
h
+
ψ
.
Let a
1
=
(
ϕ
+
η
u
2
)
/
(
μ
h
+
ψ
), b
1
=
1
-
u
1, c
=
β
v
h
ϑ
α
h, d
=
κ
1
θ
α
h, g
=
μ
h
+
α
h, h
=
β
v
h
ϑ, and k
=
μ
h
+
ϕ
+
η
u
2
+
δ
h. Then substituting r in the i equation we obtain
(30)
i
=
b
1
(
1
-
a
1
i
)
c
z
+
d
+
g
+
b
h
z
b
1
c
z
+
g
+
b
1
h
k
z
,
(
b
1
c
z
+
g
+
b
1
h
k
z
+
b
1
c
a
1
z
b
1
c
z
+
g
+
b
1
h
k
z
)
i
=
b
1
c
z
+
d
+
g
+
b
1
h
z
b
1
c
z
+
g
+
b
1
h
k
z
,
i
=
b
1
c
z
+
d
+
g
+
b
1
h
z
b
1
c
z
+
g
+
b
1
h
k
z
+
b
1
c
a
1
z
.
Using equation five of the system (24), we obtain
(31)
(
μ
v
+
α
v
+
τ
u
3
)
y
+
β
h
v
i
y
=
β
h
v
ϑ
(
1
-
z
)
i
,
y
=
(
1
-
z
)
β
h
v
ϑ
i
μ
v
+
α
v
+
τ
u
3
+
β
h
v
ϑ
i
.
Solving for z gives
(32)
α
v
y
=
(
μ
v
+
τ
u
3
)
z
,
z
=
α
v
y
μ
v
+
τ
u
3
.
Let m
2
=
α
v
β
h
v
ϑ, n
2
=
μ
v
+
τ
u
3, p
1
=
μ
v
+
α
v
+
τ
u
3, and q
1
=
β
h
v
ϑ. Substituting y in the z equation gives
(33)
z
=
(
1
-
z
)
α
v
β
h
v
ϑ
i
(
μ
v
+
τ
u
3
)
(
μ
v
+
α
v
+
τ
u
3
+
β
h
v
ϑ
i
)
=
m
2
i
n
2
(
p
1
+
q
1
i
)
-
m
2
i
z
n
2
(
p
1
+
q
1
i
)
=
m
2
i
n
2
p
1
+
n
2
q
1
i
+
m
2
i
.
Substituting z in the i equation gives
(34)
i
=
d
+
g
+
(
b
1
c
+
b
1
h
)
(
m
2
i
/
(
n
2
p
1
+
n
2
q
1
i
+
m
2
i
)
)
g
+
(
b
1
c
+
b
1
h
k
+
b
1
c
a
1
)
(
m
2
i
/
(
n
2
p
1
+
n
2
q
1
i
+
m
2
i
)
)
.
The existence of the endemic equilibrium in Φ, can be determined when i
≠
0. After some algebraic manipulations we obtain
(35)
A
i
2
+
B
i
+
C
=
0
,
where
(36)
A
=
g
n
2
q
1
+
g
m
2
+
b
1
c
m
2
+
b
1
h
k
m
2
+
b
1
c
a
1
m
2
=
(
μ
h
+
α
h
)
(
μ
v
+
τ
u
3
)
β
h
v
ϑ
+
(
μ
h
+
α
h
)
α
v
β
h
v
ϑ
+
(
1
-
u
1
)
β
v
h
β
h
v
α
h
α
v
ϑ
2
+
(
1
-
u
1
)
(
μ
h
+
ϕ
+
η
u
2
+
δ
h
)
β
v
h
β
h
v
α
v
ϑ
2
+
(
1
-
u
1
)
(
ϕ
+
η
u
2
μ
h
+
ψ
)
β
v
h
β
h
v
α
h
α
v
ϑ
2
>
0
,
C
=
-
d
n
2
p
1
-
n
2
p
1
=
-
(
μ
v
+
τ
u
3
)
(
μ
v
+
α
v
+
τ
u
3
)
[
(
1
-
u
1
)
]
<
0
,
B
=
d
n
2
q
1
+
d
m
2
+
g
n
2
q
1
+
g
m
2
+
b
1
h
m
2
+
g
n
2
p
1
-
b
1
c
m
2
=
(
μ
v
+
τ
u
3
)
κ
1
θ
α
h
β
h
v
ϑ
+
κ
1
θ
α
h
α
v
β
h
v
ϑ
+
(
μ
h
+
α
h
)
(
μ
v
+
τ
u
3
)
β
h
v
ϑ
+
(
μ
h
+
α
h
)
α
v
β
h
v
ϑ
+
(
1
-
u
1
)
β
v
h
β
h
v
α
v
ϑ
2
+
(
μ
h
+
α
h
)
(
μ
v
+
τ
u
3
)
(
μ
v
+
α
v
+
τ
u
3
)
-
β
h
v
β
v
h
ϑ
2
α
h
α
v
(
1
-
u
1
)
=
(
μ
v
+
τ
u
3
)
κ
1
θ
α
h
β
h
v
ϑ
+
κ
1
θ
α
h
α
v
β
h
v
ϑ
+
(
μ
h
+
α
h
)
(
μ
v
+
τ
u
3
)
β
h
v
ϑ
+
(
μ
h
+
α
h
)
α
v
β
h
v
ϑ
+
(
1
-
u
1
)
β
v
h
β
h
v
α
v
ϑ
2
+
(
μ
h
+
α
h
)
(
μ
v
+
τ
u
3
)
(
μ
v
+
α
v
+
τ
u
3
)
(
1
-
R
e
2
)
.
For R
e
>
1, the existence of endemic equilibria is determined by the presence of positive real solutions of the quadratic expression (35). Consider
(37)
B
=
(
μ
v
+
τ
u
3
)
κ
1
θ
α
h
β
h
v
ϑ
+
κ
1
θ
α
h
α
v
β
h
v
ϑ
+
(
μ
h
+
α
h
)
(
μ
v
+
τ
u
3
)
β
h
v
ϑ
+
(
μ
h
+
α
h
)
κ
2
θ
α
v
β
h
v
ϑ
+
(
1
-
u
1
)
β
v
h
β
h
v
α
v
ϑ
2
+
(
μ
h
+
α
h
)
(
μ
v
+
τ
u
3
)
×
(
μ
v
+
α
v
+
τ
u
3
)
(
1
-
R
e
2
)
<
0
.
Since C
<
0 and A
>
0, then C
/
A
<
0. Thus, there exists exactly one positive endemic equilibrium for i
∈
(
0,1
] whenever R
e
>
1. This gives the threshold for the endemic persistence. Therefore, we have proved the existence and uniqueness of the endemic equilibrium, E
1, for the system (1). This result is summarized in the following theorem.Theorem 3.
IfR
e
>
1, then the model system (23) has a unique endemic equilibrium E
1.The result in Theorem3 indicates the impossibility of backward bifurcation in the malaria model, since it has no endemic equilibrium when R
e
<
1. Thus, the model (1) has a globally asymptotically stable disease-free equilibrium whenever R
e
≤
1.Next we need to investigate the global stability property of the endemic equilibrium for the case when the disease does not induce death.


## 2.1.3. Global Stability of the Endemic Equilibrium forδ
h
=
0

The model system (1), when δ
h
=
0, has a unique endemic equilibrium. By letting
(38)
Φ
2
=
{
(
S
h
,
E
h
,
I
h
,
R
h
,
S
v
,
E
v
,
I
v
)
∈
R
+
7
:
E
h
=
I
h
l
l
l
l
=
E
v
=
I
v
=
0
(
S
h
,
E
h
,
I
h
,
R
h
,
S
v
,
E
v
,
I
v
)
∈
R
+
7
}
,
with
R
e
0
=
R
e
∣
δ
h
=
0
,
then the following can be claimed.Theorem 4.
The endemic equilibrium point of the malaria model (1) with δ
h
=
0 is globally asymptotically stable whenever R
e
0
>
1 in Φ and Φ
2.Proof.
As for the case of Theorem3, it can be shown that the unique endemic equilibrium for this special case exists only if R
e
0
>
1. Additionally, N
h
=
Λ
h
/
μ
h as t
→
∞. Letting S
h
=
Λ
h
/
μ
h
-
E
h
-
I
h
-
R
h and S
v
=
Λ
v
/
μ
v
-
E
v
-
I
v and substituting into (1) give the limiting system
(39)
d
E
h
d
t
=
(
1
-
u
1
)
λ
h
(
Λ
h
μ
h
-
E
h
-
I
h
-
R
h
)
+
κ
1
θ
Λ
h
μ
h
-
(
α
h
+
μ
h
)
E
h
,
d
I
h
d
t
=
α
h
E
h
-
(
ϕ
+
η
u
2
+
μ
h
+
δ
h
)
I
h
,
d
E
v
d
t
=
λ
v
(
Λ
v
μ
v
-
E
v
-
I
v
)
-
(
α
v
+
μ
v
+
τ
u
3
)
E
v
,
d
I
v
d
t
=
α
v
E
v
-
(
μ
v
+
τ
u
3
)
I
v
.
Dulac’s multiplier 1
/
E
h
I
h
E
v
I
v (see [19]) gives
(40)
∂
∂
E
h
[
(
1
-
u
1
)
β
v
h
ϑ
E
h
I
h
E
v
Λ
h
/
μ
h
(
Λ
h
μ
h
-
E
h
-
I
h
-
R
h
)
+
κ
1
θ
Λ
h
E
h
I
h
E
v
I
v
μ
h
-
(
α
h
+
μ
h
)
I
h
E
v
I
v
]
+
∂
∂
I
h
[
α
h
I
h
E
v
I
v
-
(
ϕ
+
η
u
2
+
μ
h
+
δ
h
)
E
h
E
v
I
v
]
+
∂
∂
E
v
[
(
α
v
+
μ
v
+
τ
u
3
)
E
h
I
h
E
v
β
h
v
ϑ
E
h
E
v
I
v
Λ
h
/
μ
h
(
Λ
v
μ
v
-
E
v
-
I
v
)
ddddddldd
-
(
α
v
+
μ
v
+
τ
u
3
)
E
h
I
h
E
v
]
+
∂
∂
I
v
[
α
v
E
h
I
h
I
v
-
(
μ
v
+
τ
u
3
)
E
h
I
h
E
v
]
=
-
β
v
h
ϑ
E
h
2
I
h
E
v
-
β
v
h
ϑ
μ
h
E
h
2
E
v
(
1
-
R
h
Λ
h
/
μ
h
)
-
κ
1
θ
Λ
h
E
h
2
I
h
E
v
I
v
μ
h
-
α
h
I
h
2
E
v
I
v
-
β
h
v
ϑ
E
h
E
v
2
I
v
[
1
-
μ
h
Λ
h
]
-
α
v
E
h
I
h
I
v
2
<
0
,
since E
h
+
I
h
+
R
h
<
Λ
h
/
μ
h and E
v
+
I
v
<
Λ
v
/
μ
v in Φ
2.
Thus, by Dulac’s criterion, there are no periodic orbits inΦ and Φ
2. Since Φ
2 is positively invariant and the endemic equilibrium exists whenever R
e
0
>
1, then from the Poincare-Bendixson Theorem [20], it follows that all solutions of the limiting system originating in Φ remain in Φ for all t. Furthermore, the absence of periodic orbits in Φ implies that the unique endemic equilibrium of the special case of the malaria model is globally asymptotically stable whenever R
e
0
>
1.The malaria model has a locally asymptotically stable disease-free equilibrium wheneverR
e
<
1 and a unique endemic equilibrium whenever R
e
>
1. In addition, the unique endemic equilibrium is globally asymptotically stable for the case δ
h
=
0 if R
e
0
>
1.In the next section, we apply the optimal control method using Pontryagin’s Maximum Principle to determine the necessary conditions for the combined optimal control of ITNs, IRS, and treatment effort which are being practised in Karonga District, Malawi.


## 3. Analysis of Optimal Control of the Malaria Model

The force of infection in the human population is reduced by a factor of(
1
-
u
1
(
t
)
) where u
1
(
t
) represents the use of ITNs as a means of minimizing or eliminating mosquito-human contact. u
2
(
t
)
∈
[
0
,
1
] represents the control effort on treatment of infectious individuals. This indeed represents the situation when individuals in the community seek treatment after visiting the hospitals or dispensary in their areas. For the mosquito population, we have a third control variable, u
3
(
t
). IRS affects the whole mosquito population by increasing its mortality rate by u
3
(
t
). We will use an approach similar to that in Lashari and Zaman [21] which consists of applying Pontryagin’s Maximum Principal to determine the conditions under which eradication of the disease can be achieved in finite time. Following the dynamics of the model system (1) with appropriate initial conditions, the bounded Lebesgue measurable control is used with the objective functional defined as
(41)
Γ
(
u
1
,
u
2
,
u
3
)
=
∫
0
T
f
‍
(
1
2
C
1
E
h
+
C
2
I
h
+
C
3
N
v
f
f
f
f
f
f
+
1
2
(
A
1
u
1
2
+
A
2
u
2
2
+
A
3
u
3
2
)
)
d
tsubject to the differential equations in (1), where C
1, C
2, C
3, A
1, A
2, A
3 are positive weights. We choose a quadratic cost on the controls in line with what is known in the literature on epidemic controls [21–23].The purpose of this section is to minimize the cost functional (41). This functional includes the exposed and infectious human population and the total mosquito population. In addition, it has the cost of implementing personal protection using ITNs, A
1
u
1
2, treatment of infected individuals, A
2
u
2
2, and spraying of houses, A
3
u
3
2. A linear function has been chosen for the cost incurred by exposed individuals, C
1
E
h, infected individuals, C
2
I
h, and the mosquito population, C
3
N
v. A quadratic form is used for the cost on the controls A
1
u
1
2, A
2
u
2
2, and A
3
u
3
2, such that the terms (
1
/
2
)
A
1
u
1
2, (
1
/
2
)
A
2
u
2
2, and (
1
/
2
)
A
3
u
3
2 describe the cost associated with the ITNs, treatment and IRS, respectively. We select to model the control efforts via a linear combination of quadratic terms u
i
2
(
t
), (
i
=
1,2
,
3
) and the constants C
i and A
i, (
i
=
1,2
,
3
) representing a measure of the relative cost of the interventions over the time horizon (
0
,
T
f
). We seek an optimal control u
1
*
(
t
), u
2
*
(
t
), and u
3
*
(
t
) such that
(42)
Γ
(
u
1
*
,
u
2
*
,
u
3
*
)
=
min
⁡
(
u
1
,
u
2
,
u
3
)
∈
Φ
2
Γ
(
u
1
,
u
2
,
u
3
)
,
where
(43)
Φ
2
=
{
[
0
,
T
f
]
u
=
(
u
1
,
u
2
,
u
3
)
∣
u
i
(
t
)
is  Lebesgue measurable
,
f
f
i
0
≤
u
(
t
)
≤
1
for
t
∈
[
0
,
T
f
]
⟶
[
0,1
]
,
i
=
1,2
,
3
}
is the control set, subject to the system (1) and appropriate initial conditions. We develop the optimal system for which the necessary conditions it must satisfy come from Pontryagin’s Maximum Principle [24].


### 3.1. Existence of an Optimal Control Problem

Pontryagin’s Maximum Principal converts (1) and (41) into a problem of minimizing pointwise the Lagrangian, L, and Hamiltonian, H, with respect to u
1, u
2, and u
3. The Lagrangian of the control problem is given by
(44)
L
=
C
1
E
h
+
C
2
I
h
+
C
3
N
v
+
1
2
(
A
1
u
1
2
+
A
2
u
2
2
+
A
3
u
3
2
)
.
We search for the minimal value of the Lagrangian. This can be done by defining the Hamiltonian, H, for the control problem as
(45)
H
=
L
+
λ
1
d
S
h
d
t
+
λ
2
d
E
h
d
t
+
λ
3
d
I
h
d
t
+
λ
4
d
R
h
d
t
+
λ
5
d
S
v
d
t
+
λ
6
d
E
v
d
t
+
λ
7
d
I
v
d
t
=
C
1
E
h
+
C
2
I
h
+
C
3
N
v
+
1
2
(
A
1
u
1
2
+
A
2
u
2
2
+
A
3
u
3
2
)
+
λ
1
[
Λ
h
+
(
1
-
κ
1
)
θ
N
h
+
(
ϕ
+
η
u
2
)
(
1
-
ρ
)
I
h
ffffffifffi
-
(
1
-
u
1
)
β
v
h
ϑ
I
v
S
h
-
μ
h
S
h
+
ψ
R
h
]
+
λ
2
[
(
1
-
u
1
)
β
v
h
ϑ
I
v
S
h
+
κ
1
θ
N
h
-
(
α
h
+
μ
h
)
E
h
]
+
λ
3
[
α
h
E
h
-
(
ϕ
+
η
u
2
+
μ
h
+
δ
h
)
I
h
]
+
λ
4
[
(
ϕ
+
η
u
2
)
ρ
I
h
-
(
μ
h
+
ψ
)
R
h
]
+
λ
5
[
Λ
v
-
β
h
v
ϑ
I
h
S
v
-
(
μ
v
+
τ
u
3
)
S
v
]
+
λ
6
[
β
h
v
ϑ
I
h
S
v
-
(
α
v
+
μ
v
+
τ
u
3
)
E
v
]
+
λ
7
[
α
v
E
v
-
(
μ
v
+
τ
u
3
)
I
v
]
.Next, we prove the existence of an optimal control for the model system (1).Theorem 5.
The model system (1) with the initial conditions at t
=
0 has control strategies and there exists an optimal control u
→
*
=
(
u
1
*
,
u
2
*
,
u
3
*
)
∈
Φ
2 such that
(46)
min
⁡
(
u
1
,
u
2
,
u
3
)
∈
Φ
2
Γ
(
u
1
,
u
2
,
u
3
)
=
Γ
(
u
1
*
,
u
2
*
,
u
3
*
)
.Proof.
The state and the control variables of the system (1) are nonnegative values. The control set Φ
2 is closed and convex. The integrand of the objective cost function Γ expressed by (1) is a convex function of (
u
1
,
u
2
,
u
3
) on the control set Φ
2. 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
(47)
Γ
(
u
1
,
u
2
,
u
3
)
≥
ξ
1
(
|
u
1
|
2
+
|
u
2
|
2
,
|
u
3
|
2
)
ϵ
/
2
-
ξ
2
.
This concludes the existence of an optimal control because the state variables are bounded.


### 3.2. Classification of the Optimal Control Problem

We use Pontryagin’s Maximum Principle to develop the necessary conditions for this optimal control since there exists an optimal control for minimizing the functional (41) subject to the system of equations in (1). From Lashari and Zaman [21], if (
χ
,
u
) is an optimal solution of an optimal control problem, then there exists a nontrivial vector function λ
=
(
λ
1
,
λ
2
,
λ
3
,
…
,
λ
n
) satisfying the following equations:
(48)
0
=
∂
H
(
t
,
χ
,
u
,
λ
)
∂
u
λ
′
=
∂
H
(
t
,
χ
,
u
,
λ
)
∂
χ
d
χ
d
t
=
-
∂
H
(
t
,
χ
,
u
,
λ
)
∂
λ
.
Hence the necessary conditions of the Hamiltonian, H, can be applied in (45).Theorem 6.
For the optimal control triple(
u
1
*
,
u
2
*
,
u
3
*
) with their optimal state solutions (
S
h
*
,
E
h
*
,
I
h
*
,
R
h
*
,
S
v
*
,
E
v
*
,
I
v
*
) that minimizes Γ
(
u
1
,
u
2
,
u
3
) over Φ
2, then there exist adjoint variables (
λ
1
,
λ
2
,
λ
3
,
λ
4
,
λ
5
,
λ
6
,
λ
7
) satisfying
(49)
-
λ
1
′
=
(
1
-
κ
1
)
θ
λ
1
+
(
1
-
u
1
)
β
v
h
ϑ
(
λ
2
-
λ
1
)
I
v
-
μ
h
λ
1
+
θ
κ
1
λ
2
,
-
λ
2
′
=
(
1
-
κ
1
)
θ
λ
1
+
θ
κ
1
λ
2
+
α
h
(
λ
3
-
λ
2
)
-
μ
h
λ
2
+
C
1
,
-
λ
3
′
=
(
1
-
κ
1
)
θ
λ
1
+
(
ϕ
+
η
u
2
)
[
(
1
-
ρ
)
λ
1
+
ρ
λ
4
]
+
θ
κ
1
λ
2
-
(
ϕ
+
η
u
2
+
μ
h
+
δ
h
)
λ
3
+
β
h
v
ϑ
(
λ
6
-
λ
5
)
S
v
+
C
2
,
-
λ
4
′
=
(
1
-
κ
1
)
θ
λ
1
+
ψ
λ
1
+
θ
κ
1
λ
2
-
(
μ
h
+
ψ
)
λ
4
,
-
λ
5
′
=
β
h
v
ϑ
(
λ
6
-
λ
5
)
I
h
-
(
μ
v
+
τ
u
3
)
λ
5
+
C
3
,
-
λ
6
′
=
α
v
(
λ
7
-
λ
6
)
-
(
μ
v
+
τ
u
3
)
λ
6
+
C
3
,
-
λ
7
′
=
(
1
-
u
1
)
β
v
h
ϑ
(
λ
2
-
λ
1
)
S
h
-
(
μ
v
+
τ
u
3
)
λ
7
+
C
3
with transversality conditions
(50)
λ
1
(
T
f
)
=
λ
2
(
T
f
)
=
λ
3
(
T
f
)
=
λ
4
(
T
f
)
=
λ
5
(
T
f
)
=
λ
6
(
T
f
)
=
λ
7
(
T
f
)
=
0
.
Additionally, the optimal control triple (
u
1
*
,
u
2
*
,
u
3
*
) that minimizes Γ over Φ
2 satisfies the optimality condition
(51)
u
1
*
=
max
⁡
{
0
,
min
⁡
(
1
,
β
v
h
ϑ
(
λ
2
-
λ
1
)
I
v
*
S
h
*
A
1
)
}
,
u
2
*
=
max
⁡
{
0
,
min
⁡
(
1
,
η
(
λ
3
+
ρ
(
λ
1
-
λ
4
)
-
λ
1
)
I
h
*
A
2
)
}
,
u
3
*
=
max
⁡
{
0
,
min
⁡
(
1
,
τ
(
λ
5
S
v
*
+
λ
6
E
v
*
+
λ
7
I
v
*
)
A
3
)
}
.Proof.
The adjoint equations can be determined by using the differential equations governing the adjoint variables. The Hamiltonian function,H, is differentiated with respect to S
h, E
h, I
h, R
h, S
v, E
v, and I
v. The adjoint equation is given by
(52)
-
d
λ
1
d
t
=
∂
H
∂
S
h
=
(
1
-
κ
1
)
θ
λ
1
+
(
1
-
u
1
)
β
v
h
ϑ
(
λ
2
-
λ
1
)
I
v
-
μ
h
λ
1
+
θ
κ
1
λ
2
,
-
d
λ
2
d
t
=
∂
H
∂
E
h
=
(
1
-
κ
1
)
θ
λ
1
+
θ
κ
1
λ
2
+
α
h
(
λ
3
-
λ
2
)
-
μ
h
λ
2
+
C
1
,
-
d
λ
3
d
t
=
∂
H
∂
I
h
=
(
1
-
κ
1
)
θ
λ
1
+
(
ϕ
+
η
u
2
)
[
(
1
-
ρ
)
λ
1
+
ρ
λ
4
]
+
θ
κ
1
λ
2
-
(
ϕ
+
η
u
2
+
μ
h
+
δ
h
)
λ
3
+
β
h
v
ϑ
(
λ
6
-
λ
5
)
S
v
+
C
2
,
-
d
λ
4
d
t
=
∂
H
∂
R
h
=
(
1
-
κ
1
)
θ
λ
1
+
ψ
λ
1
+
θ
κ
1
λ
2
-
(
μ
h
+
ψ
)
λ
4
,
-
d
λ
5
d
t
=
∂
H
∂
S
v
=
β
h
v
ϑ
(
λ
6
-
λ
5
)
I
h
-
(
μ
v
+
τ
u
3
)
λ
5
+
C
3
,
-
d
λ
6
d
t
=
∂
H
∂
E
v
=
α
v
(
λ
7
-
λ
6
)
-
(
μ
v
+
τ
u
3
)
λ
6
+
C
3
,
-
d
λ
7
d
t
=
∂
H
∂
I
v
=
(
1
-
u
1
)
β
v
h
ϑ
(
λ
2
-
λ
1
)
S
h
-
(
μ
v
+
τ
u
3
)
λ
7
+
C
3
,
with the transversality conditions
(53)
λ
1
(
T
f
)
=
λ
2
(
T
f
)
=
λ
3
(
T
f
)
=
λ
4
(
T
f
)
=
λ
5
(
T
f
)
=
λ
6
(
T
f
)
=
λ
7
(
T
f
)
=
0
.
Solving ∂
H
/
∂
u
1
=
0, ∂
H
/
∂
u
2
=
0, and ∂
H
/
∂
u
3
=
0, evaluating at the optimal control on the interior of the control set, where 0
<
u
i
<
1, for i
=
1,2
,
3, and letting S
h
=
S
h
*, E
h
=
E
h
*, I
h
=
I
h
*, R
h
=
R
h
*, S
v
=
S
v
*, E
v
=
E
v
*, and I
v
=
I
v
* yields
(54)
∂
H
∂
u
1
=
A
u
1
*
+
β
v
h
ϑ
λ
1
I
v
*
S
h
*
-
β
v
h
ϑ
λ
2
I
v
*
S
h
*
=
0
,
∂
H
∂
u
2
=
A
2
u
2
+
(
1
-
ρ
)
η
λ
1
I
h
*
-
η
λ
3
I
h
*
+
η
ρ
λ
4
I
h
*
=
0
,
∂
H
∂
u
3
=
A
3
u
3
*
-
τ
λ
5
S
v
*
-
τ
λ
6
E
v
*
-
τ
λ
7
I
v
*
=
0
,
for which
(55)
u
1
*
=
β
v
h
ϑ
(
λ
2
-
λ
1
)
I
v
*
S
h
*
A
1
,
u
2
*
=
η
(
λ
3
+
ρ
(
λ
1
-
λ
4
)
-
λ
1
)
I
h
*
A
2
,
u
3
=
τ
(
λ
5
S
v
*
+
λ
6
E
v
*
+
λ
7
I
v
*
)
A
3
.
Then
(56)
u
1
*
=
max
⁡
{
0
,
min
⁡
(
1
,
β
v
h
ϑ
(
λ
2
-
λ
1
)
I
v
*
S
h
*
A
1
)
}
,
u
2
*
=
max
⁡
{
0
,
min
⁡
(
1
,
η
(
λ
3
+
ρ
(
λ
1
-
λ
4
)
-
λ
1
)
I
h
*
A
2
)
}
,
u
3
*
=
max
⁡
{
0
,
min
⁡
(
1
,
τ
(
λ
5
S
v
*
+
λ
6
E
v
*
+
λ
7
I
v
*
)
A
3
)
}
.
We achieve the uniqueness of the optimal control for small T
f due to the prior boundedness of the state and adjoint functions and the resulting Lipschitz structure of the ordinary differential equations. The uniqueness of the optimal control triple trails from the uniqueness of the optimal system, which consists of (1), (49), and (50) with characterization of the optimal control (51).The optimality system is comprised of the state system (1), the adjoint system (49), initial conditions at t
=
0, boundary conditions (50), and the characterization of the optimal control (51). Hence the state and optimal control can be calculated using the optimality system. Hence using the fact that the second derivatives of the Lagrangian with respect to u
1, u
2, and u
3, respectively, are positive indicates that the optimal problem is a minimum at controls u
1
*, u
2
*, and u
3
*. Substituting u
1
*, u
2
*, and u
3
* in the system (1), we obtain
(57)
d
S
h
*
d
t
=
Λ
h
+
(
1
-
κ
1
)
θ
N
h
*
d
S
h
*
d
t
=
+
(
ϕ
+
η
u
2
)
(
1
-
ρ
)
I
h
*
-
β
v
h
ϑ
I
v
*
S
h
*
d
S
h
*
d
t
=
×
(
1
-
max
⁡
{
0
,
min
⁡
(
1
,
β
v
h
ϑ
(
λ
2
-
λ
1
)
I
v
*
S
h
*
A
1
)
}
)
d
S
h
*
d
t
=
-
μ
h
S
h
*
+
ψ
R
h
*
,
d
E
h
*
d
t
=
β
v
h
ϑ
I
v
*
S
h
*
d
E
h
*
d
t
=
×
(
1
-
max
⁡
{
0
,
min
⁡
(
1
,
β
v
h
ϑ
(
λ
2
-
λ
1
)
I
v
*
S
h
*
A
1
)
}
)
d
E
h
*
d
t
=
+
κ
1
θ
N
h
*
-
α
h
E
h
*
-
μ
h
E
h
*
,
d
I
h
*
d
t
=
α
h
E
h
*
-
(
(
max
⁡
{
0
,
min
⁡
(
1
,
η
(
λ
3
+
ρ
(
λ
1
-
λ
4
)
-
λ
1
)
I
h
*
A
2
)
}
)
ϕ
+
η
a
a
i
a
a
×
(

max
⁡
{
0
,
min
⁡
(
1
,
η
(
λ
3
+
ρ
(
λ
1
-
λ
4
)
-
λ
1
)
I
h
*
A
2
)
}

)
a
i
d
d
a
a
l
a
a
+
μ
h
+
δ
h
(

max
⁡
{
0
,
min
⁡
(
1
,
η
(
λ
3
+
ρ
(
λ
1
-
λ
4
)
-
λ
1
)
I
h
*
A
2
)
}
)
)
I
h
*
,
d
R
h
*
d
t
=
(
ϕ
+
η
(
max
⁡
{
I
h
*
×
(
A
2
)
-
1
)
0
,
min
⁡
(
(
A
2
)
-
1
1
,
η
(
λ
3
+
ρ
(
λ
1
-
λ
4
)
-
λ
1
)
a
a
a
a
a
a
a
a
a
d
d
d
d
d
d
d
d
a
a
a
a
a
×
I
h
*
×
(
A
2
)
-
1
)
}
)
)
ρ
I
h
dddddll
-
(
μ
h
+
ψ
)
R
h
,
d
S
v
*
d
t
=
Λ
v
-
(
(
A
2
)
-
1
μ
v
+
τ
ffffffdddf
×
(
max
⁡
{
0
,
min
⁡
(
×
(
A
3
)
-
1
1
,
τ
(
λ
5
S
v
*
+
λ
6
E
v
*
+
λ
7
I
v
*
)
ffifffffffffdddldddddfffffiffi
×
(
A
3
)
-
1
)
}
)
)
S
v
*
-
β
h
v
ϑ
I
h
*
S
v
*
,
d
E
v
*
d
t
=
β
h
v
ϑ
I
h
*
S
v
*
-
(
μ
v
+
τ
(
max
⁡
{
×
(
A
3
)
-
1
)
0
,
min
⁡
(
×
(
A
3
)
-
1
1
,
τ
(
λ
5
S
v
*
+
λ
6
E
v
*
+
λ
7
I
v
*
)
ffddddddddffffffddddddffffi
×
(
A
3
)
-
1
)
}
)
)
E
v
*
-
α
v
E
v
*
,
d
I
v
*
d
t
=
α
v
E
v
-
(
μ
v
+
τ
(
max
⁡
{
×
(
A
3
)
-
1
)
0
,
min
⁡
(
×
(
A
3
)
-
1
1
,
τ
(
λ
5
S
v
*
+
λ
6
E
v
*
+
λ
7
I
v
*
)
z
z
z
z
z
z
z
d
d
d
d
d
d
d
f
z
z
f
f
z
×
(
A
3
)
-
1
)
}
)
)
I
v
*
with H
* at (
t
,
S
h
*
,
E
h
*
,
I
h
*
,
R
h
*
,
u
1
*
,
u
2
*
,
u
3
*
,
λ
1
,
λ
2
,
…
,
λ
7
):
(58)
H
*
=
C
1
E
h
+
C
2
I
h
+
C
3
N
v
+
1
2
(
A
1
(
max
⁡
{
0
,
min
⁡
(
1
,
β
v
h
ϑ
(
λ
2
-
λ
1
)
I
v
*
S
h
*
A
1
)
}
)
2
fffffi
+
A
2
(
max
⁡
{
(
A
2
)
-
1
0
,
min
⁡
(
(
A
2
)
-
1
1
,
η
(
λ
3
+
ρ
(
λ
1
-
λ
4
)
-
λ
1
)
I
h
*
f
ffdddddddddfdfffffffffflff
×
(
A
2
)
-
1
)
}
)
2
fffffi

+

A
3
(
max
⁡
{
0
,
min
⁡
(
1
,
τ
(
λ
5
S
v
*
+
λ
6
E
v
*
+
λ
7
I
v
*
)
A
3
)
}

)
2

(
max
⁡
{
0
,
min
⁡
(
1
,
β
v
h
ϑ
(
λ
2
-
λ
1
)
I
v
*
S
h
*
A
1
)
}
)
2
)
+
λ
1
d
S
h
*
d
t
+
λ
2
d
E
h
*
d
t
+
λ
3
d
I
h
*
d
t
+
λ
4
d
R
h
*
d
t
+
λ
5
d
S
v
*
d
t
+
λ
6
d
E
v
*
d
t
+
λ
7
d
I
v
*
d
t
.
We solve the systems (57) and (58) numerically to determine the optimal control and the state.


### 3.3. Numerical Results on Optimal Control Analysis

In this section we discuss the method and present the results obtained from solving the optimality system numerically using the parameter values in Table11.Here we consider the optimal control values of the three intervention strategies, namely, ITNs, IRS, and treatment, which are common strategies in Karonga District, Malawi. In Figure2(a) we see that if the three intervention strategies are effectively implemented and used, they have a positive impact compared to having ITNs and IRS (u
1 and u
3), respectively, as the only intervention strategies in the community. The initial increase in the infected human population in the graph with interventions of ITNs and IRS may be due to the fact that some people refused to have their houses sprayed with insecticide chemicals due to their primitive traditional beliefs. In addition, as Karonga District is along the shore of Lake Malawi, some members of the community do not use ITNs owing to negative beliefs in the chemicals used and also due to hot weather in the districts. On the other hand, the district is waterlogged and hence this leads to an increase in mosquitoes breeding sites.Simulations of model (1) showing the effects of intervention measures.
(a)
(b)
(c)
(d)Similar results appear in Figure2(b) where the impact of the intervention strategies are compared with the use of ITNs (u
1) and treatment (u
2). The results show that the concurrent administered intervention strategies lead to a decrease in the number of infected human population much faster than when ITNs and treatment are used as the only intervention strategies in the community. A similar occurrence is observed when IRS (u
3) and treatment (u
2) are used as the only means of intervention strategies in the community (see Figure 2(c)).In addition, we also looked at the effects of these intervention measures as a stand-alone approach of preventing or controlling malaria disease in the community. Figure2(d) depicts the comparison of the effects of each intervention strategy and the effect of a multi-intervention strategy. The figure indicates that if the three combined intervention measures are effectively practised, the infected human population is much lower compared to a situation where only one intervention strategy is used. The graph of treatment (u
2) practised as the only means of intervention measure shows a high number of infected human population owing to a number of reasons. One of the reasons is that, in this situation, the mosquito population is unaffected; hence the infected mosquitoes will still be available in the community causing more infections to susceptible humans. Furthermore, most people in the area do not visit the dispensary or hospital for medication when they observe signs or symptoms of malaria disease since they need to cover long distances to reach the hospital. The interviews conducted revealed that some individuals opt to simply use the medication left by the previous patient or they buy medicines from the shops and use it before being diagnosed. Consequently treatment needs to be consolidated with preventive measures such as ITNs and IRS for optimal control.The epidemiological implication of the above result is that malaria could be eliminated from the community if prevention and treatment can lead to a situation whereR
e is less than unity. However, other factors need to be considered.


## 3.1. Existence of an Optimal Control Problem

Pontryagin’s Maximum Principal converts (1) and (41) into a problem of minimizing pointwise the Lagrangian, L, and Hamiltonian, H, with respect to u
1, u
2, and u
3. The Lagrangian of the control problem is given by
(44)
L
=
C
1
E
h
+
C
2
I
h
+
C
3
N
v
+
1
2
(
A
1
u
1
2
+
A
2
u
2
2
+
A
3
u
3
2
)
.
We search for the minimal value of the Lagrangian. This can be done by defining the Hamiltonian, H, for the control problem as
(45)
H
=
L
+
λ
1
d
S
h
d
t
+
λ
2
d
E
h
d
t
+
λ
3
d
I
h
d
t
+
λ
4
d
R
h
d
t
+
λ
5
d
S
v
d
t
+
λ
6
d
E
v
d
t
+
λ
7
d
I
v
d
t
=
C
1
E
h
+
C
2
I
h
+
C
3
N
v
+
1
2
(
A
1
u
1
2
+
A
2
u
2
2
+
A
3
u
3
2
)
+
λ
1
[
Λ
h
+
(
1
-
κ
1
)
θ
N
h
+
(
ϕ
+
η
u
2
)
(
1
-
ρ
)
I
h
ffffffifffi
-
(
1
-
u
1
)
β
v
h
ϑ
I
v
S
h
-
μ
h
S
h
+
ψ
R
h
]
+
λ
2
[
(
1
-
u
1
)
β
v
h
ϑ
I
v
S
h
+
κ
1
θ
N
h
-
(
α
h
+
μ
h
)
E
h
]
+
λ
3
[
α
h
E
h
-
(
ϕ
+
η
u
2
+
μ
h
+
δ
h
)
I
h
]
+
λ
4
[
(
ϕ
+
η
u
2
)
ρ
I
h
-
(
μ
h
+
ψ
)
R
h
]
+
λ
5
[
Λ
v
-
β
h
v
ϑ
I
h
S
v
-
(
μ
v
+
τ
u
3
)
S
v
]
+
λ
6
[
β
h
v
ϑ
I
h
S
v
-
(
α
v
+
μ
v
+
τ
u
3
)
E
v
]
+
λ
7
[
α
v
E
v
-
(
μ
v
+
τ
u
3
)
I
v
]
.Next, we prove the existence of an optimal control for the model system (1).Theorem 5.
The model system (1) with the initial conditions at t
=
0 has control strategies and there exists an optimal control u
→
*
=
(
u
1
*
,
u
2
*
,
u
3
*
)
∈
Φ
2 such that
(46)
min
⁡
(
u
1
,
u
2
,
u
3
)
∈
Φ
2
Γ
(
u
1
,
u
2
,
u
3
)
=
Γ
(
u
1
*
,
u
2
*
,
u
3
*
)
.Proof.
The state and the control variables of the system (1) are nonnegative values. The control set Φ
2 is closed and convex. The integrand of the objective cost function Γ expressed by (1) is a convex function of (
u
1
,
u
2
,
u
3
) on the control set Φ
2. 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
(47)
Γ
(
u
1
,
u
2
,
u
3
)
≥
ξ
1
(
|
u
1
|
2
+
|
u
2
|
2
,
|
u
3
|
2
)
ϵ
/
2
-
ξ
2
.
This concludes the existence of an optimal control because the state variables are bounded.


## 3.2. Classification of the Optimal Control Problem

We use Pontryagin’s Maximum Principle to develop the necessary conditions for this optimal control since there exists an optimal control for minimizing the functional (41) subject to the system of equations in (1). From Lashari and Zaman [21], if (
χ
,
u
) is an optimal solution of an optimal control problem, then there exists a nontrivial vector function λ
=
(
λ
1
,
λ
2
,
λ
3
,
…
,
λ
n
) satisfying the following equations:
(48)
0
=
∂
H
(
t
,
χ
,
u
,
λ
)
∂
u
λ
′
=
∂
H
(
t
,
χ
,
u
,
λ
)
∂
χ
d
χ
d
t
=
-
∂
H
(
t
,
χ
,
u
,
λ
)
∂
λ
.
Hence the necessary conditions of the Hamiltonian, H, can be applied in (45).Theorem 6.
For the optimal control triple(
u
1
*
,
u
2
*
,
u
3
*
) with their optimal state solutions (
S
h
*
,
E
h
*
,
I
h
*
,
R
h
*
,
S
v
*
,
E
v
*
,
I
v
*
) that minimizes Γ
(
u
1
,
u
2
,
u
3
) over Φ
2, then there exist adjoint variables (
λ
1
,
λ
2
,
λ
3
,
λ
4
,
λ
5
,
λ
6
,
λ
7
) satisfying
(49)
-
λ
1
′
=
(
1
-
κ
1
)
θ
λ
1
+
(
1
-
u
1
)
β
v
h
ϑ
(
λ
2
-
λ
1
)
I
v
-
μ
h
λ
1
+
θ
κ
1
λ
2
,
-
λ
2
′
=
(
1
-
κ
1
)
θ
λ
1
+
θ
κ
1
λ
2
+
α
h
(
λ
3
-
λ
2
)
-
μ
h
λ
2
+
C
1
,
-
λ
3
′
=
(
1
-
κ
1
)
θ
λ
1
+
(
ϕ
+
η
u
2
)
[
(
1
-
ρ
)
λ
1
+
ρ
λ
4
]
+
θ
κ
1
λ
2
-
(
ϕ
+
η
u
2
+
μ
h
+
δ
h
)
λ
3
+
β
h
v
ϑ
(
λ
6
-
λ
5
)
S
v
+
C
2
,
-
λ
4
′
=
(
1
-
κ
1
)
θ
λ
1
+
ψ
λ
1
+
θ
κ
1
λ
2
-
(
μ
h
+
ψ
)
λ
4
,
-
λ
5
′
=
β
h
v
ϑ
(
λ
6
-
λ
5
)
I
h
-
(
μ
v
+
τ
u
3
)
λ
5
+
C
3
,
-
λ
6
′
=
α
v
(
λ
7
-
λ
6
)
-
(
μ
v
+
τ
u
3
)
λ
6
+
C
3
,
-
λ
7
′
=
(
1
-
u
1
)
β
v
h
ϑ
(
λ
2
-
λ
1
)
S
h
-
(
μ
v
+
τ
u
3
)
λ
7
+
C
3
with transversality conditions
(50)
λ
1
(
T
f
)
=
λ
2
(
T
f
)
=
λ
3
(
T
f
)
=
λ
4
(
T
f
)
=
λ
5
(
T
f
)
=
λ
6
(
T
f
)
=
λ
7
(
T
f
)
=
0
.
Additionally, the optimal control triple (
u
1
*
,
u
2
*
,
u
3
*
) that minimizes Γ over Φ
2 satisfies the optimality condition
(51)
u
1
*
=
max
⁡
{
0
,
min
⁡
(
1
,
β
v
h
ϑ
(
λ
2
-
λ
1
)
I
v
*
S
h
*
A
1
)
}
,
u
2
*
=
max
⁡
{
0
,
min
⁡
(
1
,
η
(
λ
3
+
ρ
(
λ
1
-
λ
4
)
-
λ
1
)
I
h
*
A
2
)
}
,
u
3
*
=
max
⁡
{
0
,
min
⁡
(
1
,
τ
(
λ
5
S
v
*
+
λ
6
E
v
*
+
λ
7
I
v
*
)
A
3
)
}
.Proof.
The adjoint equations can be determined by using the differential equations governing the adjoint variables. The Hamiltonian function,H, is differentiated with respect to S
h, E
h, I
h, R
h, S
v, E
v, and I
v. The adjoint equation is given by
(52)
-
d
λ
1
d
t
=
∂
H
∂
S
h
=
(
1
-
κ
1
)
θ
λ
1
+
(
1
-
u
1
)
β
v
h
ϑ
(
λ
2
-
λ
1
)
I
v
-
μ
h
λ
1
+
θ
κ
1
λ
2
,
-
d
λ
2
d
t
=
∂
H
∂
E
h
=
(
1
-
κ
1
)
θ
λ
1
+
θ
κ
1
λ
2
+
α
h
(
λ
3
-
λ
2
)
-
μ
h
λ
2
+
C
1
,
-
d
λ
3
d
t
=
∂
H
∂
I
h
=
(
1
-
κ
1
)
θ
λ
1
+
(
ϕ
+
η
u
2
)
[
(
1
-
ρ
)
λ
1
+
ρ
λ
4
]
+
θ
κ
1
λ
2
-
(
ϕ
+
η
u
2
+
μ
h
+
δ
h
)
λ
3
+
β
h
v
ϑ
(
λ
6
-
λ
5
)
S
v
+
C
2
,
-
d
λ
4
d
t
=
∂
H
∂
R
h
=
(
1
-
κ
1
)
θ
λ
1
+
ψ
λ
1
+
θ
κ
1
λ
2
-
(
μ
h
+
ψ
)
λ
4
,
-
d
λ
5
d
t
=
∂
H
∂
S
v
=
β
h
v
ϑ
(
λ
6
-
λ
5
)
I
h
-
(
μ
v
+
τ
u
3
)
λ
5
+
C
3
,
-
d
λ
6
d
t
=
∂
H
∂
E
v
=
α
v
(
λ
7
-
λ
6
)
-
(
μ
v
+
τ
u
3
)
λ
6
+
C
3
,
-
d
λ
7
d
t
=
∂
H
∂
I
v
=
(
1
-
u
1
)
β
v
h
ϑ
(
λ
2
-
λ
1
)
S
h
-
(
μ
v
+
τ
u
3
)
λ
7
+
C
3
,
with the transversality conditions
(53)
λ
1
(
T
f
)
=
λ
2
(
T
f
)
=
λ
3
(
T
f
)
=
λ
4
(
T
f
)
=
λ
5
(
T
f
)
=
λ
6
(
T
f
)
=
λ
7
(
T
f
)
=
0
.
Solving ∂
H
/
∂
u
1
=
0, ∂
H
/
∂
u
2
=
0, and ∂
H
/
∂
u
3
=
0, evaluating at the optimal control on the interior of the control set, where 0
<
u
i
<
1, for i
=
1,2
,
3, and letting S
h
=
S
h
*, E
h
=
E
h
*, I
h
=
I
h
*, R
h
=
R
h
*, S
v
=
S
v
*, E
v
=
E
v
*, and I
v
=
I
v
* yields
(54)
∂
H
∂
u
1
=
A
u
1
*
+
β
v
h
ϑ
λ
1
I
v
*
S
h
*
-
β
v
h
ϑ
λ
2
I
v
*
S
h
*
=
0
,
∂
H
∂
u
2
=
A
2
u
2
+
(
1
-
ρ
)
η
λ
1
I
h
*
-
η
λ
3
I
h
*
+
η
ρ
λ
4
I
h
*
=
0
,
∂
H
∂
u
3
=
A
3
u
3
*
-
τ
λ
5
S
v
*
-
τ
λ
6
E
v
*
-
τ
λ
7
I
v
*
=
0
,
for which
(55)
u
1
*
=
β
v
h
ϑ
(
λ
2
-
λ
1
)
I
v
*
S
h
*
A
1
,
u
2
*
=
η
(
λ
3
+
ρ
(
λ
1
-
λ
4
)
-
λ
1
)
I
h
*
A
2
,
u
3
=
τ
(
λ
5
S
v
*
+
λ
6
E
v
*
+
λ
7
I
v
*
)
A
3
.
Then
(56)
u
1
*
=
max
⁡
{
0
,
min
⁡
(
1
,
β
v
h
ϑ
(
λ
2
-
λ
1
)
I
v
*
S
h
*
A
1
)
}
,
u
2
*
=
max
⁡
{
0
,
min
⁡
(
1
,
η
(
λ
3
+
ρ
(
λ
1
-
λ
4
)
-
λ
1
)
I
h
*
A
2
)
}
,
u
3
*
=
max
⁡
{
0
,
min
⁡
(
1
,
τ
(
λ
5
S
v
*
+
λ
6
E
v
*
+
λ
7
I
v
*
)
A
3
)
}
.
We achieve the uniqueness of the optimal control for small T
f due to the prior boundedness of the state and adjoint functions and the resulting Lipschitz structure of the ordinary differential equations. The uniqueness of the optimal control triple trails from the uniqueness of the optimal system, which consists of (1), (49), and (50) with characterization of the optimal control (51).The optimality system is comprised of the state system (1), the adjoint system (49), initial conditions at t
=
0, boundary conditions (50), and the characterization of the optimal control (51). Hence the state and optimal control can be calculated using the optimality system. Hence using the fact that the second derivatives of the Lagrangian with respect to u
1, u
2, and u
3, respectively, are positive indicates that the optimal problem is a minimum at controls u
1
*, u
2
*, and u
3
*. Substituting u
1
*, u
2
*, and u
3
* in the system (1), we obtain
(57)
d
S
h
*
d
t
=
Λ
h
+
(
1
-
κ
1
)
θ
N
h
*
d
S
h
*
d
t
=
+
(
ϕ
+
η
u
2
)
(
1
-
ρ
)
I
h
*
-
β
v
h
ϑ
I
v
*
S
h
*
d
S
h
*
d
t
=
×
(
1
-
max
⁡
{
0
,
min
⁡
(
1
,
β
v
h
ϑ
(
λ
2
-
λ
1
)
I
v
*
S
h
*
A
1
)
}
)
d
S
h
*
d
t
=
-
μ
h
S
h
*
+
ψ
R
h
*
,
d
E
h
*
d
t
=
β
v
h
ϑ
I
v
*
S
h
*
d
E
h
*
d
t
=
×
(
1
-
max
⁡
{
0
,
min
⁡
(
1
,
β
v
h
ϑ
(
λ
2
-
λ
1
)
I
v
*
S
h
*
A
1
)
}
)
d
E
h
*
d
t
=
+
κ
1
θ
N
h
*
-
α
h
E
h
*
-
μ
h
E
h
*
,
d
I
h
*
d
t
=
α
h
E
h
*
-
(
(
max
⁡
{
0
,
min
⁡
(
1
,
η
(
λ
3
+
ρ
(
λ
1
-
λ
4
)
-
λ
1
)
I
h
*
A
2
)
}
)
ϕ
+
η
a
a
i
a
a
×
(

max
⁡
{
0
,
min
⁡
(
1
,
η
(
λ
3
+
ρ
(
λ
1
-
λ
4
)
-
λ
1
)
I
h
*
A
2
)
}

)
a
i
d
d
a
a
l
a
a
+
μ
h
+
δ
h
(

max
⁡
{
0
,
min
⁡
(
1
,
η
(
λ
3
+
ρ
(
λ
1
-
λ
4
)
-
λ
1
)
I
h
*
A
2
)
}
)
)
I
h
*
,
d
R
h
*
d
t
=
(
ϕ
+
η
(
max
⁡
{
I
h
*
×
(
A
2
)
-
1
)
0
,
min
⁡
(
(
A
2
)
-
1
1
,
η
(
λ
3
+
ρ
(
λ
1
-
λ
4
)
-
λ
1
)
a
a
a
a
a
a
a
a
a
d
d
d
d
d
d
d
d
a
a
a
a
a
×
I
h
*
×
(
A
2
)
-
1
)
}
)
)
ρ
I
h
dddddll
-
(
μ
h
+
ψ
)
R
h
,
d
S
v
*
d
t
=
Λ
v
-
(
(
A
2
)
-
1
μ
v
+
τ
ffffffdddf
×
(
max
⁡
{
0
,
min
⁡
(
×
(
A
3
)
-
1
1
,
τ
(
λ
5
S
v
*
+
λ
6
E
v
*
+
λ
7
I
v
*
)
ffifffffffffdddldddddfffffiffi
×
(
A
3
)
-
1
)
}
)
)
S
v
*
-
β
h
v
ϑ
I
h
*
S
v
*
,
d
E
v
*
d
t
=
β
h
v
ϑ
I
h
*
S
v
*
-
(
μ
v
+
τ
(
max
⁡
{
×
(
A
3
)
-
1
)
0
,
min
⁡
(
×
(
A
3
)
-
1
1
,
τ
(
λ
5
S
v
*
+
λ
6
E
v
*
+
λ
7
I
v
*
)
ffddddddddffffffddddddffffi
×
(
A
3
)
-
1
)
}
)
)
E
v
*
-
α
v
E
v
*
,
d
I
v
*
d
t
=
α
v
E
v
-
(
μ
v
+
τ
(
max
⁡
{
×
(
A
3
)
-
1
)
0
,
min
⁡
(
×
(
A
3
)
-
1
1
,
τ
(
λ
5
S
v
*
+
λ
6
E
v
*
+
λ
7
I
v
*
)
z
z
z
z
z
z
z
d
d
d
d
d
d
d
f
z
z
f
f
z
×
(
A
3
)
-
1
)
}
)
)
I
v
*
with H
* at (
t
,
S
h
*
,
E
h
*
,
I
h
*
,
R
h
*
,
u
1
*
,
u
2
*
,
u
3
*
,
λ
1
,
λ
2
,
…
,
λ
7
):
(58)
H
*
=
C
1
E
h
+
C
2
I
h
+
C
3
N
v
+
1
2
(
A
1
(
max
⁡
{
0
,
min
⁡
(
1
,
β
v
h
ϑ
(
λ
2
-
λ
1
)
I
v
*
S
h
*
A
1
)
}
)
2
fffffi
+
A
2
(
max
⁡
{
(
A
2
)
-
1
0
,
min
⁡
(
(
A
2
)
-
1
1
,
η
(
λ
3
+
ρ
(
λ
1
-
λ
4
)
-
λ
1
)
I
h
*
f
ffdddddddddfdfffffffffflff
×
(
A
2
)
-
1
)
}
)
2
fffffi

+

A
3
(
max
⁡
{
0
,
min
⁡
(
1
,
τ
(
λ
5
S
v
*
+
λ
6
E
v
*
+
λ
7
I
v
*
)
A
3
)
}

)
2

(
max
⁡
{
0
,
min
⁡
(
1
,
β
v
h
ϑ
(
λ
2
-
λ
1
)
I
v
*
S
h
*
A
1
)
}
)
2
)
+
λ
1
d
S
h
*
d
t
+
λ
2
d
E
h
*
d
t
+
λ
3
d
I
h
*
d
t
+
λ
4
d
R
h
*
d
t
+
λ
5
d
S
v
*
d
t
+
λ
6
d
E
v
*
d
t
+
λ
7
d
I
v
*
d
t
.
We solve the systems (57) and (58) numerically to determine the optimal control and the state.


## 3.3. Numerical Results on Optimal Control Analysis

In this section we discuss the method and present the results obtained from solving the optimality system numerically using the parameter values in Table11.Here we consider the optimal control values of the three intervention strategies, namely, ITNs, IRS, and treatment, which are common strategies in Karonga District, Malawi. In Figure2(a) we see that if the three intervention strategies are effectively implemented and used, they have a positive impact compared to having ITNs and IRS (u
1 and u
3), respectively, as the only intervention strategies in the community. The initial increase in the infected human population in the graph with interventions of ITNs and IRS may be due to the fact that some people refused to have their houses sprayed with insecticide chemicals due to their primitive traditional beliefs. In addition, as Karonga District is along the shore of Lake Malawi, some members of the community do not use ITNs owing to negative beliefs in the chemicals used and also due to hot weather in the districts. On the other hand, the district is waterlogged and hence this leads to an increase in mosquitoes breeding sites.Simulations of model (1) showing the effects of intervention measures.
(a)
(b)
(c)
(d)Similar results appear in Figure2(b) where the impact of the intervention strategies are compared with the use of ITNs (u
1) and treatment (u
2). The results show that the concurrent administered intervention strategies lead to a decrease in the number of infected human population much faster than when ITNs and treatment are used as the only intervention strategies in the community. A similar occurrence is observed when IRS (u
3) and treatment (u
2) are used as the only means of intervention strategies in the community (see Figure 2(c)).In addition, we also looked at the effects of these intervention measures as a stand-alone approach of preventing or controlling malaria disease in the community. Figure2(d) depicts the comparison of the effects of each intervention strategy and the effect of a multi-intervention strategy. The figure indicates that if the three combined intervention measures are effectively practised, the infected human population is much lower compared to a situation where only one intervention strategy is used. The graph of treatment (u
2) practised as the only means of intervention measure shows a high number of infected human population owing to a number of reasons. One of the reasons is that, in this situation, the mosquito population is unaffected; hence the infected mosquitoes will still be available in the community causing more infections to susceptible humans. Furthermore, most people in the area do not visit the dispensary or hospital for medication when they observe signs or symptoms of malaria disease since they need to cover long distances to reach the hospital. The interviews conducted revealed that some individuals opt to simply use the medication left by the previous patient or they buy medicines from the shops and use it before being diagnosed. Consequently treatment needs to be consolidated with preventive measures such as ITNs and IRS for optimal control.The epidemiological implication of the above result is that malaria could be eliminated from the community if prevention and treatment can lead to a situation whereR
e is less than unity. However, other factors need to be considered.


## 4. Analysis of Malaria Data from Karonga District

A structured questionnaire was developed and administered in Karonga District, Malawi, in order to determine how intervention strategies of malaria disease are being practised and their effectiveness. The questionnaire was used to conduct a directed one to one interview, on the respondents who were randomly sampled, with total size of 502, which means that 502 questionnaires were administered. The enumerators received training before they went to the field, where the questionnaire was discussed with the respondents.We consider statistical results of how intervention strategies are practised. Different graphs and tables are depicted for all the prevention and treatment strategies.Knowledge of malaria transmission is a key to malaria prevention [25]. Figure 3 shows that 409 respondents (81.5%) stated that malaria transmission occurs when bitten by an infected mosquito, contrary to popular belief in developing countries (UNICEF (2000) [26]) that someone may be infected with malaria when soaked in water, where most respondents (96.8%) answered “no.” Table 4 shows the number of individuals who mentioned that malaria disease and malaria transmission are through a bite from an infected mosquito. 79.9% responded “yes” to both malaria being a major disease and malaria transmission being acquired through a bite from an infected mosquito.Table 4
Major diseases in the community-malaria (malaria transmission-bite from an infected mosquito cross tabulation).
Malaria transmission-bite from an infected mosquito
Total
No
Yes
Major diseases in the community-malaria
No
Count
11
15
26
% of total
2.2%
3.0%
5.3%
Yes
Count
73
393
466
% of total
14.8%
79.9%
94%
Total
Count
84
408
492
% of total
17.1%
82.9%
100.0%Figure 3
Knowledge of malaria transmission.According to the Center for Disease Control and Prevention [27], prevention is better than cure and there are a number of methods which people can use to prevent malaria. Most of the respondents (90.4%) stated that they use ITNs or long lasting insecticide treated bed nets (LLITNs) as a method of preventing occurrence of malaria. Only 16.5% of the respondents mentioned that their houses are sprayed (IRS) (see Figure 4). From Table 5, there are about 1% cases of malaria occurrences and about 2.2% respondents had not experienced any malaria cases after spraying their houses with IRS. Some of the respondents used nets as well as spraying their houses. Table 6 shows that 35.5% of the respondents had experienced malaria occurrence after their houses were sprayed and also used the bed nets. From the respondents who had been using ITNs or LLITNs, 44.9% did not have any malaria occurrence. However, about 45.9% had had an occurrence of malaria despite the use of bed nets (see Table 7).Table 5
Malaria occurrence after spraying and without using net.
Frequency
Percent
No
11
2.2
Yes
5
1.0
No response
486
96.8
Total
502
100.0Table 6
Malaria occurrence after spraying and using net.
Frequency
Percent
No
195
38.8
Yes
178
35.5
No response
129
25.7
Total
502
100.0Table 7
Method of preventing malaria use of ITN or LLITN (use of net and malaria occurrence cross tabulation).
Use of net and malaria
Total
Occurrence
No
Yes
Method of preventing use of ITN or LLITN
No
Count
17
27
44
% of total
3.6%
5.7%
9.2%
Yes
Count
214
219
433
% of total
44.9%
45.9%
90.8%
Total
Count
231
246
477
% of total
48.4%
51.6%
100.0%Figure 4
Methods of preventing malaria.A chi-square test of independence was conducted with a chi-square value of 1.454 since the cross tabulation is a two by two table. The assumption of no cells having an expected value less than 5 was not violated, hence the use of chi-square test. Table8 shows that a calculated value of 1.454 was obtained. The level of significance used for the chi-square test was 0.05 which was compared to a P value of 0.228. Since 0.228
>
0.05, the null hypothesis was rejected; hence, there was an association between malaria occurrence and preventive methods (use of ITNs or LLITNs) and their proportions are not significantly different. In Table 10, a indicates the number of cells with expected count less than 5 while b shows Yates continuity correction which is calculated for 2 by 2 table. The relationship of indoor residual spraying (IRS) against the occurrence of malaria after the house was sprayed and use of bed nets was examined. From Table 9 for the respondents who had not had their houses sprayed, 42.5% did not have an occurrence of malaria, while for those respondents who had their houses sprayed, 38.4% had an occurrence of malaria. A chi-square test of independence was conducted with a chi-square value of 0.019. The assumption of no cells having an expected value less than 5 was not violated, hence the use of chi-square test. The level of significance used for the chi-square test was 0.05 which was compared to a P value of 0.889. Table 10 shows a P value of 0.889
>
0.05, which means there was an association between preventive method (IRS) and malaria occurrence after spraying and using ITNs but, however, there were no significant differences in proportions between preventive methods and malaria occurrence after spraying.Table 8
Chi-square test of use of ITN or LLITN and use of net and malaria occurrence.
Value
df
Asymp. sig. (2-sided)
Pearson chi-square
1.86
1
a
1
0.173
Continuity correctionb
1.454
1
0.228Table 9
Method of preventing malaria house sprayed IRS (malaria occurrence after spraying and using net cross tabulation).
Malaria occurrence after spraying and using net
Total
No
Yes
Method of preventing house spread (IRS)
No
Count
158
143
301
% of total
42.5%
38.4%
80.9%
Yes
Count
36
35
71
% of total
9.7%
9.4%
19.1%
Total
Count
194
178
372
% of total
52.2%
47.8%
100.0%Table 10
House sprayed IRS: malaria occurrence after spraying and using net.
Value
df
Asymp. sig. (2-sided)
Pearson chi-square
0.074
a
1
0.786
Continuity correctionb
0.019
1
0.889Table 11
Table for parameter values of the optimal malaria model.
Symbol
Value
Source
Λ
h
0.02326
Calculated
θ
0.3
Calculated
μ
v
0.1429
[6]
α
h
1/17
[7]
β
h
v
0.0375
Assumed
ρ
0.035
Assumed
Ψ
0.0018
Assumed
u
1
(
t
)
0.0904
Calculated
u
3
(
t
)
0.076
Calculated
η
0.4
Assumed
N
h
20,000
[8]
κ
1
0.003
Calculated
μ
h
0.04326
[8]
δ
h
0.03454
[8]
β
v
h
0.0833
[9]
ϑ
0.35
Assumed
ϕ
0.3
Assumed
α
v
0.1
Assumed
u
2
(
t
)
0.165
Calculated
τ
0.01
Assumed
Λ
v
1000 day−1
[7]
N
v
3,500
AssumedDespite knowledge of malaria transmission and prevention, malaria cases still occur [28]. About 50% of the respondents who used ITNs mentioned that malaria still occurs. This is probably because bed nets are only used when going to bed, hence the vulnerability. Furthermore the proportion of those respondents who used ITNs or LLITNs and suffered from malaria was not significantly different from those who did not use ITNs or LLITNs but suffer from malaria (χ
2
=
1.454, calculated P value = 0.228, level of significance of 0.05). Of the respondents who had had their houses sprayed (about 3.2%), 1% experienced an occurrence of malaria. Those respondents who had their houses sprayed and suffered from malaria even though they had used ITNs or LLITNs showed no significant difference with those respondents who had had their houses sprayed and had not suffered from malaria even though they had used ITNs or LLITNs (χ
2
=
0.019, calculated P value = 0.889, level of significance of 0.05).


## 5. Numerical Results

The numerical simulations and analysis were carried out using a fourth order Runge-Kutta scheme in Matlab. Our aim was to determine and verify the analytic results and the stability of the model system (1). Some of the parameter values were calculated from the data collected in Karonga District, Malawi, between the months of January and September, 2013. The other parameter values were obtained from the National Statistical Office (NSO) in Zomba, Malawi, some have been assumed, and very few have been taken from the literature. The Government of Malawi has organized a number of intervention strategies in order to fight against malaria in the country through the National Malaria Control Program (NMCP) [29].


### 5.1. Dynamics of Human State Variables for the Malaria Model without Intervention Strategies

The analysis of the model without intervention strategies was carried out in order to determine the dynamics of the disease in the population. The simulation was generated in a four-year time frame since the first campaign of malaria intervention strategies in Karonga District was performed in the year 2010. The susceptible human population is decreasing exponentially (see Figure5(a)) showing that most susceptible humans are exposed to the disease due to unavailability of intervention strategies. This has led to an exponential increase in the exposed human population (Figure 5(b)) and the infected population (Figure 5(c)) with R
0
=
1.8894. The infected human population increases due to an increase in the exposure of susceptible individuals toPlasmodium falciparum. This means that thePlasmodium falciparum will continue to multiply in the human and mosquito populations since there are no intervention strategies to reduce or eradicate the disease. Hence, there is a need of having intervention strategies in order to reduce or eradicate the disease.The changes in the four state variables of the malaria model without intervention strategies illustrating the dynamics, with time, of (a) susceptible individuals, (b) exposed individuals, (c) infected individuals, and (d) recovered human individuals.
(a)
(b)
(c)
(d)


### 5.2. Prevalence in the Malaria Model without Intervention Strategies

Prevalence is defined as the ratio of the number of cases of the disease in a population to the total number of individuals in population at a given time. The disease prevalence in Figure6 shows a steady increase during the first days of infection due to a high number of individuals withoutPlasmodium falciparum. The graph drops asymptotically due to a reduced number of susceptible human population, showing evidence that communities can be wiped out with malaria disease if none of the intervention methods are speedily put into place.Figure 6
Exposed humans without any intervention strategies.


### 5.3. Dynamical System of the Individual Population State Variables in the Model with Intervention Strategies

We now consider the effects of the three intervention strategies (ITNs, IRS, and treatment) which are campaigned concurrently in Karonga District, Malawi. It appears that if the three intervention strategies are effectively monitored and implemented, then there is a positive impact of combating malaria in the community. In Figures7(a) and 7(b), an increase is observed in the human population recovering compared with the exposed and infected individuals. This steady increase in the recovery class is due to the readily available effective prevention strategies and treatment. Figure 7(a), where human population is plotted against time on a four-year period, shows a steady increase in the recovery of the human population. It is evidenced that during this period chemicals which are available in the mosquito nets and sprayed in the houses are effective in reducing the mosquito population and at the same time reducing the contact rate between the human population and the mosquito population. Treatment strategy has also played an important role in reducing the number of infected individuals thus leading to an increase in recovered individuals.The dynamics of the exposed and infected individuals in the model with intervention strategies forR
e
=
0.6922.
(a)
(b)A slight decrease in the gradient of the recovered human population is observed in Figure7(b) as the period of using prevention and treatment strategies available is increased. Despite being given ITNs and spraying long lasting chemicals in the houses, the campaign needs to be revisited to confirm whether the chemicals are still effective. The decrease of the graph explains that the effectiveness of the treated bed nets and the residual spraying deteriorates with time. Hence there is need to assess the right time interval to carry out the respraying of the houses and the resupply of ITNs.


## 5.1. Dynamics of Human State Variables for the Malaria Model without Intervention Strategies

The analysis of the model without intervention strategies was carried out in order to determine the dynamics of the disease in the population. The simulation was generated in a four-year time frame since the first campaign of malaria intervention strategies in Karonga District was performed in the year 2010. The susceptible human population is decreasing exponentially (see Figure5(a)) showing that most susceptible humans are exposed to the disease due to unavailability of intervention strategies. This has led to an exponential increase in the exposed human population (Figure 5(b)) and the infected population (Figure 5(c)) with R
0
=
1.8894. The infected human population increases due to an increase in the exposure of susceptible individuals toPlasmodium falciparum. This means that thePlasmodium falciparum will continue to multiply in the human and mosquito populations since there are no intervention strategies to reduce or eradicate the disease. Hence, there is a need of having intervention strategies in order to reduce or eradicate the disease.The changes in the four state variables of the malaria model without intervention strategies illustrating the dynamics, with time, of (a) susceptible individuals, (b) exposed individuals, (c) infected individuals, and (d) recovered human individuals.
(a)
(b)
(c)
(d)


## 5.2. Prevalence in the Malaria Model without Intervention Strategies

Prevalence is defined as the ratio of the number of cases of the disease in a population to the total number of individuals in population at a given time. The disease prevalence in Figure6 shows a steady increase during the first days of infection due to a high number of individuals withoutPlasmodium falciparum. The graph drops asymptotically due to a reduced number of susceptible human population, showing evidence that communities can be wiped out with malaria disease if none of the intervention methods are speedily put into place.Figure 6
Exposed humans without any intervention strategies.


## 5.3. Dynamical System of the Individual Population State Variables in the Model with Intervention Strategies

We now consider the effects of the three intervention strategies (ITNs, IRS, and treatment) which are campaigned concurrently in Karonga District, Malawi. It appears that if the three intervention strategies are effectively monitored and implemented, then there is a positive impact of combating malaria in the community. In Figures7(a) and 7(b), an increase is observed in the human population recovering compared with the exposed and infected individuals. This steady increase in the recovery class is due to the readily available effective prevention strategies and treatment. Figure 7(a), where human population is plotted against time on a four-year period, shows a steady increase in the recovery of the human population. It is evidenced that during this period chemicals which are available in the mosquito nets and sprayed in the houses are effective in reducing the mosquito population and at the same time reducing the contact rate between the human population and the mosquito population. Treatment strategy has also played an important role in reducing the number of infected individuals thus leading to an increase in recovered individuals.The dynamics of the exposed and infected individuals in the model with intervention strategies forR
e
=
0.6922.
(a)
(b)A slight decrease in the gradient of the recovered human population is observed in Figure7(b) as the period of using prevention and treatment strategies available is increased. Despite being given ITNs and spraying long lasting chemicals in the houses, the campaign needs to be revisited to confirm whether the chemicals are still effective. The decrease of the graph explains that the effectiveness of the treated bed nets and the residual spraying deteriorates with time. Hence there is need to assess the right time interval to carry out the respraying of the houses and the resupply of ITNs.


## 6. Conclusions

An optimal control model (using a deterministic system of nonlinear ordinary differential equations) for the transmission dynamics of malaria in Karonga District, Malawi, was presented. The model considered a varying total human population that incorporated recruitment of new individuals into the susceptible class through birth or immigration, and those immigrant individuals who were exposed to the disease were recruited into the exposed individual class. The prevention (IRS and ITNs) and other treatment intervention strategies were included in the model to assess the potential impact of these strategies on the transmission dynamics of the disease.Our model incorporated features that were potentially effective to control or reduce the transmission of malaria disease in Malawi. Analysis of the optimal control model revealed that there exists a domain where the model is epidemiologically and mathematically well-posed. We also computed the effective reproduction number,R
e, and then qualitatively analyzed the existence and stability of their model equilibria. The basic reproduction number, R
0, was obtained from the threshold reproduction number by eliminating all the intervention strategies. Then it was proved that if R
e
<
1, the disease cannot survive in the district. Hence the effective reproduction number, R
e, is an essential indication of the effort required to eliminate the disease. It was also found that R
e
≤
R
0 which implied that increased preventive and control intervention practices had a positive impact on the reduction of R
e. Thus, malaria can be eradicated in the district by deployment of a combination of intervention strategies such as effective mass drug administration and vector control (LLITNs and IRS) to combat and eventually eliminate the disease.Analysis of the model supported that effective control or eradication of malaria can be achieved by the combination of protection and treatment measures. We have seen that when the three intervention strategies are combined, there is a greater reduction in the number of exposed and infected individuals. The prevention strategies played a greater role in reducing the number of infected individuals by lowering the contact rate between the mosquito and human populations, for instance, through the use of ITNs. On the other hand both prevention strategies led to the reduction of the mosquito population hence lowering the infected mosquito population. Effective treatment consolidated the prevention strategies. This study provides useful tools for assessing the effectiveness of a combination of the three intervention strategies and analyzing the potential impact of prevention with treatment.Demographic findings also showed that the preventive measures, namely, ITNs and IRS, if effectively practised, can help to reduce malaria transmission. These primary health intervention strategies are very important as they reduce the mosquito population and contact between the human and mosquito populations. These practices will lead to a reduction in the transfer ofPlasmodium between the host and the vector. However, Dzinjalamala [30] states that malaria control in Malawi is still heavily reliant on chemotherapy. Hence the approach needs to change and effectively accommodate the campaign strategy taking place in Karonga District in order to combat the disease. Therefore, the presumptive treatment for fever and the primary health intervention practices (LLITNs and IRS) should both be effectively implemented or practised in order to reduce or eliminate malaria disease.Using the optimal values of the three intervention practices (LLITNs, IRS, and treatment), the results showed that the combination of the three intervention strategies has a positive and greater impact in eliminating or reducing the epidemic of malaria. This can be achieved when the measures are effectively implemented by the suppliers and effectively practised by the beneficiaries (the community members).To effectively control and potentially eradicate the spread of malaria, treatment programs must be complemented with other intervention strategies such as vector reduction and personal protection. Intervention practices that involve both prevention and treatment controls yield relatively better results. The combination of these strategies can play a positive role in Karonga District in reducing or eradicating malaria disease. Therefore, control and prevention efforts aimed at lowering the infectivity of infected individuals to the mosquito vector will contribute greatly to the reduction of malaria transmission and this will eventually lower the prevalence of malaria and the incidence of the disease in the community.The proposed model has some limitations. We did not consider infective immigrants. Also, the population was not stratified by age as it is well-known that malaria disproportionately affects children under the age of 5 years. Finally, there are various control measures out there and only three basic ones were considered herein.



---

*Source: 594256-2014-06-19.xml*