Mathematical assessment of the role of vector insecticide resistance and feeding/resting behavior on malaria transmission dynamics: Optimal control analysis.

The large-scale use of insecticide-treated bednets (ITNs) and indoor residual spraying (IRS), over the last two decades, has resulted in a dramatic reduction of malaria incidence globally. However, the effectiveness of these interventions is now being threatened by numerous factors, such as resistance to insecticide in the mosquito vector and their preference to feed and rest outdoors or early in the evening (when humans are not protected by the bednets). This study presents a new deterministic model for assessing the population-level impact of mosquito insecticide resistance on malaria transmission dynamics. A notable feature of the model is that it stratifies the mosquito population in terms of type (wild or resistant to insecticides) and feeding preference (indoor or outdoor). The model is rigorously analysed to gain insight into the existence and asymptotic stability properties of the various disease-free equilibria of the model namely the trivial disease-free equilibrium, the non-trivial resistant-only boundary disease-free equilibrium and a non-trivial disease-free equlibrium where both the wild and resistant mosquito geneotypes co-exist). Simulations of the model, using data relevant to malaria transmission dynamics in Ethiopia (a malaria-endemic nation), show that the use of optimal ITNs alone, or in combination with optimal IRS, is more effective than the singular implementation of an optimal IRS-only strategy. Further, when the effect of the fitness cost of insecticide resistance with respect to fecundity (i.e., assuming a decrease in the baseline birth rate of new resistant-type adult female mosquitoes) is accounted for, numerical simulations of the model show that the combined optimal ITNs-IRS strategy could lead to the effective control of the disease, and insecticide resistance effectively managed during the first 8 years of the 15-year implementation period of the insecticides-based anti-malaria control measures in the community.

Introduction

Malaria, caused by Plasmodium parasites, is a major life-threatening disease (WHO, 2016). The parasites are transmitted to humans through the bites of infected adult female Anopheles mosquitoes. In the year 2015, for example, 91 countries and regions in sub-Saharan Africa, South-East Asia, Latin America and the Middle East had ongoing malaria transmission and nearly half of the world's population was at risk of malaria (WHO, 2016). The people at highest risk of malaria infection include infants, children under 5 years of age, pregnant women and people living with HIV/AIDS (owing to their weak or not fully-developed immune system) (Mohammed-Awel & Numfor, 2017; WHO, 2016). Over two-thirds (70%) of all malaria deaths occur in children under age of five (WHO). According to the 2016 report of the World Health Organization (WHO), there were 212 million cases of malaria in 2015 and 429,000 deaths. In places where the mosquito lifespan is longer (such as some places in Africa (WHO, 2016)), transmission is higher. This is because the parasite has time to complete its development cycle (sporogonic) inside the mosquito (WHO). The lifespan of the African Anopheles mosquito species is long, and it has strong human-biting habit (due to these and other reasons, approximately 90% of the cases, and 92% of the deaths, occurred in the African region in 2015 (WHO, 2016)). Malaria transmission also depends on climatic conditions, such as rainfall patterns, temperature and humidity (Okuneye & Gumel, 2017; WHO, 2016; WHO) (these conditions are known to significantly affect the survival and population size (abundance) of mosquitoes).

Numerous control strategies, such as the use of insecticide-treated mosquito bednets (ITNs) or long-lasting insecticidal nets (LLINs) and indoor residual spraying (IRS) with chemical insecticides, antimalarial drugs (such as Artemisinin-based combination therapies (ACTs) (WHO, 2016)), and intermittent preventive treatment of pregnant women and infants, are being used to reduce (or prevent) malaria transmission in endemic areas. Furthermore, several candidate anti-malaria vaccines are being developed (Asale et al., 2014; Churcher, Lissenden, Griffin, Worrall, & Ranson, 2016; Mohammed-Awel, Zhao, Numfor, & Lenhart, 2017; Mohammed-Awel & Numfor, 2017; Ngonghala, Del Valle, Zhao, & Mohammed-Awel, 2014; Ngonghala, Mohammed, Zhao, & Prosper, 2016; Okumu & Moore, 2011; WHO, 2016). Owing to these efforts (i.e., implementation of ITNs and/or IRS strategies), malaria incidence and malaria mortality have been significantly reduced during the period 2010–2015 (WHO, 2016). It is now generally believed that ITNs and IRS are the most effective control strategies for reducing (or preventing) malaria transmission in endemic settings (Birget & Koella, 2015; Brown, Dickinson, & Kramer, 2013; Chitnis, Hyman, & Cushing, 2008; Choi et al., 2014; Churcher et al., 2016; Corbel et al., 2010; Griffin et al., 2010; Jones et al., 2013; Ngonghala et al., 2014; WHO, 2016). LLINs are a form of ITNs specifically designed to remain effective for many years without re-treatment (Anuse, Sahu, Subramanian, & Gunasekaran, 2015). The use of ITNs or IRS is promoted as the major control strategy in the malaria control and elimination plan endorsed by WHO (Churcher et al., 2016; WHO, 2016). The two vector control methods (ITNs and IRS) are mostly used concurrently, within the same households. Some studies have suggested that ITNs or IRS, if used singly, may not be enough to effectively combat malaria transmission (particularly in holoendemic and hyperendemic areas), and that these two control measures should preferably be combined in those areas (Griffin et al., 2010; Okumu & Moore, 2011; WHO, 2016). However, it remains unclear whether the concurrent use of ITNs and IRS in the same household reduces transmission significantly, relative to using either ITNs or IRS alone (WHO, 2016) (this forms one of the main modeling objectives of the current study).

Of the four major classes of chemical insecticides currently used in malaria control efforts (namely, pyrethroids, organochlorines, organophosphates and carbamates), only the pyrethroids are approved for use in LLINs (owing to their low mammalian toxicity and irritant effect on mosquitoes), while all four are used in IRS (WHO, 2016). Although the intensive use of ITNs (especially) and IRS, singly or in combination, has led to a significant decline in malaria-related morbidity in endemic areas (studies estimate 81% of the reduction of malaria burden recorded in the past 15 years are due to the use of LLINs and IRS, with LLINs accounting for most of the reductions), this widespread and heavy use of insecticides has, unfortunately, resulted in the emergence of vector resistance to nearly every currently available agent (Alout, Roche, Dabir, & Cohuet, 2017; Birget & Koella, 2015; Brown et al., 2013; Choi et al., 2014; Corbel et al., 2010; Dondorp et al., 2009; Gourley, Liu, & Wu, 2011; Griffin et al., 2010; Jones et al., 2013; Okumu & Moore, 2011; Ranson & Lissenden, 2016; White, Griffin, et al., 2011; WHO, 2016). This problem of insecticide resistance, which is more pronounced in malaria-endemic regions in sub-Saharan Africa (Choi et al., 2014; Corbel et al., 2010; Gourley et al., 2011; Griffin et al., 2010; Jones et al., 2013; Ranson & Lissenden, 2016; White, Griffin, et al., 2011; WHO, 2016), is, in fact, worsening particularly in sub-Saharan Africa (Ranson & Lissenden, 2016). Rotational use of different classes of insecticides for IRS is recommended as one approach to manage insecticide resistance (WHO, 2016). If left unchecked, the problem of insecticide resistance could lead to substantial increases in malaria incidence and mortality, with devastating public health consequences (Ranson & Lissenden, 2016; WHO, 2016). Another important feature of mosquito-borne diseases is the feeding and resting nature (and preferences) of adult female mosquitoes. While some adult female mosquitoes feed indoors (endophagic), others feed outdoors (exophagic). Furthermore, some mosquitoes prefer to rest indoors (endophilic) or outdoors (exophilic). Consequently, endophilism, facilitates the application of ITNs and IRS to control the population of adult female mosquitoes (and, therefore, control disease spread) (Gimnig et al., 1452).

Mathematical models have been widely used to assess the population-level impact and effectiveness of various control strategies against malaria (see, for instance, (Agusto et al., 2013; Anderson & May, 1991; Chitnis et al., 2008; Mohammed-Awel et al., 2017; Mohammed-Awel & Numfor, 2017; Ngonghala et al., 2014; Ngonghala et al., 2016; Prosper, Ruktanonchai, & Martcheva, 2014; White, Griffin, et al., 2011)). However, only few of these models, which incorporate ITNs and/or IRS, considered the effect of insecticide resistance (Barbosa & Hastings, 2012; Birget & Koella, 2015; Brown et al., 2013; Gourley et al., 2011; Wairimu & Ronoh, 2016). Barbosa et al. (Barbosa & Hastings, 2012) developed a genetic model to predict changes in mosquito fitness and resistance allele frequency (parameters that describe insecticide selection, fitness cost as well as ITNs and synergist (Piperolyn butoxide (PBO)) are incorporated). The results of their investigation show that resistance was most sensitive to selection coefficients, fitness cost and dominance coefficients. Birget et al. (Birget & Koella, 2015) developed a population-genetic model of the spread of insecticide-resistance in Anopheles mosquitoes in response to ITNs and larvacides. Their study shows that the use of indoor ITNs leads to less selection pressure than the use of insecticides as larvacides. Brown et al. (Brown et al., 2013) developed a mathematical model to investigate economically optimal strategies for mosquito control in the presence of insecticide resistance. Consistent with previous studies, their results show that fitness costs are the key elements in the computation of economically optimal resistance management strategies. Gourley et al. (Gourley et al., 2011) developed a mathematical model, where the adult mosquito population is split into vulnerable or resistant (based on whether the insecticide can have an effect or not). Their theoretical study, which does not incorporate fitness costs due to insecticide resistance, gives global asymptotic stability results for the non-trivial resistant-only boundary equilibrium of the model. Wairimu et al. (Wairimu & Ronoh, 2016) gave theoretical results for a mathematical model for malaria transmission with two classes of mosquitoes (sensitive and resistant to chemical insecticides) when ITNs and IRS are used as control strategies. Global asymptotic stability of the disease-free and endemic equilibria are proved.

In the current study, a new deterministic model is designed and used to assess the population-level impact of vector insecticide resistance and feeding and resting preferences on malaria transmission dynamics. The model will be used to evaluate the community-wide impact of ITNs and IRS, implemented singly or in combination. The model is formulated and fitted (using data relevant to malaria transmission dynamics in Ethiopia) in Section 2. The asymptotic stability properties of the associated disease-free (trivial and boundary) equilibria of the model are explored in 3. Optimal control analysis, based on the two controls (ITNs and IRS) is carried out in Section 4. Discussion and concluding remarks are reported in Section 5.

Model formulation

The model to be developed is for the transmission dynamics of malaria in an endemic setting which implements a control strategy based on using IRS or ITNs or their combination. It is assumed, for simplicity, that these strategies are only implemented indoors (Okumu & Moore, 2011) and vectors can travel between indoors and outdoors. The total human population at time t, denoted by , is split into the mutually-exclusive compartments of susceptible (), exposed (), infectious (), and recovered () humans, so that

The total adult female Anopheles mosquito population at time t, denoted by , is split into the total outdoor () and indoor () mosquitoes, where the population of outdoor mosquitoes is further stratified in terms of susceptible outdoor wild-type () and resistant-type (resistant to insecticides) (), exposed wild-type () and resistant-type () and infectious wild-type () and resistant-type () outdoor mosquitoes, so that

Similarly, the total population of indoor mosquitoes is classified in terms of susceptible wild-type () and resistant-type (), exposed wild-type () and resistant-type () and infectious wild-type () and resistant-type () indoor mosquitoes, so that

Hence,

Furthermore, the total wild-type outdoor (), resistant-type outdoor (), wild-type indoor (), and resistant-type indoor () mosquitoes are given, respectively, by

When ITNs are used in the community, the average number of bites per indoor mosquito per unit time (or mosquito-human contact rate indoors), denoted by , is defined as (Agusto et al., 2013; Mohammed-Awel & Numfor, 2017):where is maximum mosquito-biting rate, is minimum mosquito-biting rate, and b is insecticide-treated bednets coverage (or proportion of ITNs usage in the community).

The associated forces of infection for malaria transmission are defined by (where is the human-to-vector infection rate outdoors, is the human-to-vector infection rate indoors, is the vector-to-human infection rate outdoors, and is the vector-to-human infection rate indoors):where is the average number of bites per mosquito in outdoor per unit time, is mosquito biting rate, is transmission probability from infectious mosquitoes to susceptible humans, and is transmission probability from infectious humans to susceptible mosquitoes. It is assumed that mosquitoes are not produced indoors (i.e., it is assumed there is no suitable breeding habitats indoors for mosquitoes to breed). For the wild-type and resistant-type adult female mosquitoes outdoors, the following Verhulst-Pearl logistic growth functions, (for outdoor wild-type adult female mosquitoes) and (for outdoor resistant-type adult female mosquitoes), are chosen (Mohammed-Awel et al., 2017):where and are the production (birth) rates of new adult wild-type and resistant-type mosquitoes, respectively. Furthermore, is the environmental carrying capacity of adult female mosquitoes, and it is assumed that for all . It is further assumed that both the resistant and wild alleles are inherited (i.e., adult female mosquitoes of resistant (wild) genotype produce offsprings with the resistant (wild) genotype).

The model for the transmission dynamics of malaria, in the presence of vector control strategy based on using ITNs and IRS (and taking into account the feeding and resting preference of adult female mosquitoes (endophagic and endopholic vs. exophagic and exopholic mosquitoes)) is given by the following deterministic system of non-linear differential equations (where a dot represents differentiation with respect to time t):

A schematic diagram of the model is depicted in Fig. 1 (and the state variables and parameters of the model are described in Table 1, Table 2, respectively).

Fig. 1

The Schematic diagram of the model (2.4).

Table 1

Description of state variables of the model.

State Variable	Description
S_h	Population of susceptible humans
E_h	Population of exposed (infected but not symptomatic) humans
I_h	Population of infectious (symptomatic) humans
R_h	Population of recovered humans
S_vwo	Population of susceptible wild-type adult female outdoor mosquitoes
E_vwo	Population of exposed wild-type adult female outdoor mosquitoes
I_vwo	Population of infectious wild-type adult female outdoor mosquitoes
S_vwi	Population of susceptible wild-type adult female indoor mosquitoes
E_vwi	Population of exposed wild-type adult female indoor mosquitoes
I_vwi	Population of infectious wild-type adult female indoor mosquitoes
S_vro	Population of susceptible resistant-type adult female outdoor mosquitoes
E_vro	Population of exposed resistant-type adult female outdoor mosquitoes
I_vro	Population of infectious resistant-type adult female outdoor mosquitoes
S_vri	Population of susceptible resistant-type adult female indoor mosquitoes
E_vri	Population of exposed resistant-type adult female indoor mosquitoes
I_vri	Population of infectious resistant-type adult female indoor mosquitoes

Table 2

Description of parameters.

Parameters	Description	Baseline Value	Source
Λ_h	Human recruitment rate (due to birth or immigration)	2.45×10^6 (per year)	Estimated from CIA (2017)
μ_h	Natural death rate for humans	1/62.5 (per year)	Estimated from CIA (2017)
ν_h	Rate at which exposed humans become infectious	(1/14)×365 (per year)	Ngonghala et al. (2014)
γ_h	Recovery rate of humans	(3.5×10^−3)×365 (per year)	Ngonghala et al. (2016)
ρ_h	Rate of loss of natural immunity	(5.6×10^−3)×365 (per year)	Ngonghala et al. (2016)
δ_h	Disease-induced death rate for humans	(9.0×10^−5)×365 (per year)	Ngonghala et al. (2016)
K_v	Environmental carrying capacity of mosquitoes	4.8×10^5 (dimensionless)	Fitted
r_vw	Production (birth) rates of new adult wild-type female mosquitoes	77.4 (per year)	Fitted
r_vr	Production (birth) rates of new	76.3 (per year)	Fitted
	adult resistant-type female mosquitoes
μ_v	Natural death rate of mosquitoes	(1/14)×365 (per year)	Ngonghala et al. (2016)
σ_vw	Rate at which exposed wild-type adult resistant-type become infectious	36.5 (per year)	Ngonghala et al. 2016)
σ_vr	Rate at which exposed resistant-type mosquitoes become infectious	43.6 (per year)	Fitted
δ_iw	Death rate of wild-type mosquitoes (those exposed to insecticide) due to the use of IRS and ITNs	40.6 (per year)	Fitted
δ_ir	Death rate of resistant-type mosquitoes (those exposed to insecticide) due to the use of IRS and ITNs	15.9 (per year)	Fitted
u_i	Proportion of houses (indoors) sprayed with IRS	0.29 (dimensionless)	Estimated from National malaria program (2014)
κ	Rate of development of resistance due to the use of ITNs or IRS	10^−4 (per year)	Fitted
m_io	Mobility rate of mosquitoes from indoors to outdoors	84.99 (per year)	Fitted
m_oi	Mobility rate of mosquitoes from outdoors to indoors	77.3 (per year)	Fitted
b_max	Maximum mosquito biting rate	231.35 (per year)	Fitted
b_min	Minimum mosquito biting rate	1.0×10^−2 (per year)	Fitted
b_ho	Contact rate of mosquitoes with humans outdoors	71.44 (per year)	Fitted
β_h	Transmission probability from infectious mosquitoes to susceptible humans	22/10^3 (dimensionless)	Ngonghala et al. (2014)
β_v	Transmission probability from infectious humans to susceptible mosquitoes	48/10^2 (dimensionless)	Ngonghala et al. (2014)
b	Insecticide-treated bednets (ITNs) coverage (or proportion of ITNs usage)	0.49 (dimensionless)	Estimated from National malaria program (2014)
α_1	Modification parameter for the assumed reduction of the mobility of infectious vectors in relation to susceptible vectors (0<α_1<1)	0.9 (dimensionless)	Fitted

In the model (2.4), is the human recruitment rate (due to immigration and birth), the is rate of loss of temporary immunity acquired from prior infection (or natural immunity). Susceptible humans acquire infection following effective bites by an infected adult female Anopheles mosquito indoors (at a rate ) or outdoors (at a rate ). Humans in all epidemiological compartments are assumed to suffer natural death at a rate . Exposed humans develop clinical symptoms of malaria (and become infectious) at a rate of . Furthermore, infectious humans suffer additional death due to malaria at a rate of . Humans recover from clinical malaria at a rate .

As stated earlier, it is assumed (for simplicity) that no suitable mosquito habitats exist indoors, and that mosquito production is limited to outdoors only. The Verhulst-Pearl logistic birth function, defined in equation (2.3), is chosen for both the susceptible wild-type () and resistant-type () outdoor mosquitoes. It is assumed that resistance is inherited (that is, a resistant female adult mosquito vertically produces resistant offsprings (Gourley et al., 2011)). Outdoor susceptible mosquitoes ( and ) acquire malaria infection at the rate , and indoor susceptible mosquitoes ( and ) acquire infection at the rate . Susceptible outdoor mosquitoes move indoors at a rate , and susceptible indoor mosquitoes move outdoors at a rate . It is assumed that all mosquitoes suffer natural death at a rate of . Exposed wild-type mosquitoes ( and ) become infectious at a rate and exposed resistant-type mosquitoes ( and ) move to corresponding infectious class at a rate . Exposed outdoor mosquitoes ( and ) move indoors at a rate of , and exposed indoor mosquitoes ( and ) move outdoors at a rate of . Similarly, infectious indoor mosquitoes ( and ) move outdoors at a rate of . Infectious outdoor mosquitoes ( and ) move indoors at a rate , where the parameter accounts for the assumption that infectious mosquitoes transit (move) at a slower rate than susceptible mosquitoes. Similarly, infectious indoor mosquitoes ( and ) move outdoors at a rate .

It is further assumed that, due to the intensive use of ITNs and IRS, indoor wild-type mosquitoes develop resistance to the chemical insecticide at a rate , where κ rate of mosquito develop resistance to insecticide, is a proportion of houses (indoors) sprayed with IRS, and b is ITNs coverage (or proportion of individuals who use ITNs (Matow et al., 2015)). Female indoor mosquitoes could die when they come in contact with treated bednets or become exposed to IRS (Matow et al., 2015; Okumu & Moore, 2011). For this reason, it is assumed that indoor wild-type mosquitoes suffer additional death due to exposure to insecticides at a rate , where is death rate of wild type mosquitoes due to exposure to IRS and ITNs. Similarly, due to the use of IRS and ITNs indoors, resistant-type mosquitoes suffer additional mortality at a rate of , where is death rate of wild type mosquitoes due to the use of IRS and ITNs.

The model (2.4) is an extension of numerous malaria transmission models that assess the impact of insecticide resistance and dynamics of resistant vectors when chemical insecticides (such as ITNs, IRS, or larvacides) are used to control the vector population (such as those in (Barbosa & Hastings, 2012; Birget & Koella, 2015; Blayneh & Mohammed-Awel, 2014; Brown et al., 2013; Gourley et al., 2011; Wairimu & Ronoh, 2016)) by, inter alia,

including the use of ITNs and IRS control strategies where the vector population is stratified according to type (i.e., wild or resistant to insecticides). These classifications are not included in the genetic models in (Barbosa & Hastings, 2012; Birget & Koella, 2015; Brown et al., 2013);

incorporating mosquito feeding and resting behavior (these are not included in the models in (Barbosa & Hastings, 2012; Birget & Koella, 2015; Blayneh & Mohammed-Awel, 2014; Brown et al., 2013; Gourley et al., 2011; Wairimu & Ronoh, 2016));

carrying out optimal control analysis of the model with respect to the two insecticide-based controls (this is not done in (Barbosa & Hastings, 2012; Birget & Koella, 2015; Blayneh & Mohammed-Awel, 2014; Brown et al., 2013; Gourley et al., 2011; Wairimu & Ronoh, 2016)).

Basic properties

It is convenient to define , and The following basic properties can be established for the model (the results are standard, hence their proofs are omitted (Mohammed-Awel et al., 2017; Safi and Garba, 2012)).

All solutions of the model (2.4) with non-negative initial conditions remain non-negative for all time .

Consider the model (2.4) with . The following feasible regionis positively-invariant and attracts all solutions of the model (2.4) in .

It should be observed that the upper bound of , , is positive if (that is, if mosquito growth rate is higher than its death rate). In a closed environment, if the death rate of adult mosquitoes is higher than their birth rate, then the mosquito population eventually become extinct. In other words, the total mosquito population dies out when . Since the extinction of the mosquito population in a malaria-endemic community is ecologically unrealistic, it is assumed, from now on, that .

Data fitting and parameter estimation

The model (2.4) is fitted using yearly malaria case data for Ethiopia, for the period 2000 to 2015, extracted from (Deribew et al., 2017) (Table 3). In particular, Pearson's Chi-squared and least square method (implemented in the statistical software R) were used to fit the model (2.4) to the data. Some data parameters are obtained or estimated from literature (as described in Table 2). The model (2.4) is fitted to the data to estimate 13 unknown parameters (notably those related to insecticide resistance), namely , , , , , , κ, , , , , , and . Fig. 2 depicts the result of the data fitting (and the set of estimated/fitted parameters that best fits the model is tabulated in Table 4). To theoretically measure the goodness of the fit, the associated average relative error of the fitting was computed using the formula , where and are the exact and estimated number of cases in year (depicted in Table 3), respectively. This confirms the reasonably good fit obtained.

Table 3

Total number of new malaria cases per 100,000 (both male and female in Ethiopia) between 2000 and 2015 (extracted from (Deribew et al., 2017)).

Year	Number of new malaria cases per 100,000	Year	Number of new malaria cases per 100,000
2000	8,870	2008	4,770
2001	9,130	2009	3,050
2002	9,430	2010	1,400
2003	9,580	2011	1,300
2004	9,920	2012	1,300
2005	10,000	2013	1,130
2006	8,430	2014	1,240
2007	6,520	2015	1,180

Fig. 2

Data fitting of the model (2.4) using malaria case data from Ethiopia for the period 2000 to 2015 (given in Table 3) (Deribew et al., 2017).

Table 4

Values of estimated parameters.

Parameter	Estimated value (per year)	Parameter	Estimated value (per year)
K_v	4.8×10^8   (dimensionless)	m_io	84.99
r_vw	77.4	m_io	77.3
r_vr	76.3	b_max	231.35
σ_vr	43.6	b_min	10^−2
δ_iw	40.6	b_ho	71.44
δ_ir	15.9	α_1	0.9 (dimensionless)
κ	10^−4

Mathematical analysis

Existence of disease-free equilibria

It is convenient to, first of all, define the quantities: , , , , , , , , , , , , , , , , and .

The model (2.4) has three disease-free equilibria, namely a trivial disease-free equilibrium (TDFE; denoted by ), a non-trivial resistant-only disease-free boundary equilibrium (NTRDFE; denoted by ) and a non-trivial co-existence disease-free equilibrium (NTCDFE; denoted by ). The expressions for the three disease-free equilibria are given below:where,with,

TDFE:

NTRDFE:

NTCDFE:

Since and , it follows that and . It follows from (3.1) that NTCDFE () exists if and only if and . Similarly, it follows from (3.1) that the NTRDFE () exists only if . Furthermore,if (which is true since the reproduction fitness of the insecticide resistant mosquito reduced due to insecticide resistance).

The term in (3.3) can be expressed in the following geometric series,where (so that, is the average duration a susceptible wild-type and resistant-type mosquito spend outdoors), is the proportion of indoor mosquitoes in class that survived and moved outdoors, and is the proportion of outdoor mosquitoes in class that survived and moved indoors. Similarly, can be expressed as , where is the proportion of indoor mosquitoes in class that survived and moved outdoors, and is the proportion of outdoor mosquitoes in class that survived and moved to indoors.

Local asymptotic stability of NTRDFE () and NTCDFE ()

The trivial disease-free equilibrium (TDFE) is not ecologically realistic (since it entails having no mosquitoes in the population). Consequently, the asymptotic stability of this equilibrium is omitted. To show the local asymptotic stability of the NTRDFE (), it is convenient to define the quantity (which can be obtained by the routine application of the next generation operator method on the model (2.4) around the NTRDFE (Diekmann, Heesterbeek, & Metz, 1990; Driessche & Watmough, 2002)):where,with,

The result below follows from Theorem 2 in (Driessche & Watmough, 2002).

Consider the model (2.4) with . The resistant-only disease-free boundary equilibrium () is locally-asymptotically stable (LAS) if , and unstable if .

Similarly, the following reproduction threshold can be obtained by applying the next generation method on the model (2.4) around the co-existence disease-free equilibrium (NTCDFE):where,

Since , and , it follows that , , and . Hence, using Theorem 2 in (Driessche & Watmough, 2002), the following result is established.

Consider the model (2.4) with and . The non-trivial co-existence disease-free equilibrium () is LAS if , and unstable if .

The threshold quantity () (known as the reproduction number of the model) represents the average number of new malaria infections in humans or vectors generated by an infectious human or vector in a human or vector population where insecticide-based interventions (ITNs and IRS) are implemented (Anderson & May, 1991; Diekmann et al., 1990; Driessche & Watmough, 2002; Hethcote, 2000). The epidemiological significance of Theorem 3.2 is that the disease can be effectively controlled (when ) if the initial sizes of the sub-populations of the model are in the basin of attraction of the NTCDFE ().

Optimal control: formulation, analysis and simulations

Formulation

As stated earlier, this study assumes that IRS and ITNs are implemented indoors only (hence, these controls are not implemented outdoors). In order to determine the optimal strategy for IRS and ITNs coverage, the indoor-residual spraying rate () is re-defined as a function of time (i.e., ). Furthermore, the constant bednets coverage (b) is re-defined as a function of time (i.e., ), so that indoor the biting rate () is now re-defined as:Similarly, the forces of infection for malaria transmission are re-defined as:(while the expressions for the outdoor forces of infection, and , are as defined in Section 2). The indoor mosquito control interventions are implemented on a time interval , where T is the number of years for which the control program is implemented. The model (2.4), with the aforementioned non-constant time-dependent mosquito controls, is now given by:subject to the initial conditions:The optimal control problem seeks to minimize the number of humans infected with malaria, the number of infectious mosquitoes, the number of resistant-type mosquitoes, and the cost of implementation of ITNs and IRS controls in the community. In order to do this, the following objective functional is formulated:where is the total number of resistant-type mosquitoes, is the total number of indoor mosquitoes, and is the time interval over which the ITNs and IRS programs are implemented.

In (4.4), and (for ; ; and ) are positive constants that balance the relative importance of terms in the objective functional J. The integrand in gives the weighted number of humans infected with malaria and the total weighted number of infected mosquitoes; and the integrand in is the weighted number of resistant-type mosquitoes over the time period T. The term represents the total number of humans protected by ITNs, while the term is the total number of indoor mosquitoes affected by the ITNs and IRS. The term gives the total cost of implementing ITNs and IRS, the constants is the cost associated with the use of ITNs to protect humans from mosquito bite and is the cost associated with the use of the insecticide in ITNs and IRS to kill adult mosquitoes over the time period T.

Consequently, the optimal control problem is to find the pair , such thatwhere the set of all admissible controls () is given by:

The following result can be established using the approach in (Mohammed-Awel et al., 2017; Mohammed-Awel & Numfor, 2017).

Given the controls , there exists a positive, bounded solution to the initial value problem (4.2) and initial conditions (4.3).

Analysis: existence, characterization and uniqueness of optimal control

There exist optimal controls and that minimize the objective functional J subject to the system (4.2).

The Proof of Theorem 4.2, based on using the techniques in (Mohammed-Awel et al., 2017; Mohammed-Awel & Numfor, 2017), is given in Appendix A.

Pontryagin's Maximum principle (Pontryagin, Boltyanskii, Gamkrelidze, & Mishchenko, 2002) is used to derive necessary conditions that an optimal control must satisfy. This principle converts the problem of minimizing the objective functional (4.4) subject to the state system (4.2) into a problem of minimizing a Hamiltonian (), defined below, with respect to the controls :

The differential equations governing the adjoint variables () are obtained by differentiating partially with respect to each state variable of the system (4.2). This gives:

The resulting adjoint system is given by:subject to the transversality conditions (Pontryagin et al., 2002):

In characterizing the optimal control, the Hamiltonian () in Equation (4.5) is differentiated partially with respect to the controls b and , to get:on the interior of the control set, where

The result below follows from the Pontryagin's Maximum principle (Pontryagin et al., 2002).where,where is the maximum ITNs coverage in the community and the maximum IRS coverage in the community. Furthermore, the following result holds.

Given an optimal, , and solutions, , , , , , , , , , , , , , , , and of the corresponding state system (4.2), then there exist adjoint variables for , which satisfy the adjoint system in (4.6) and transversality conditions (4.7). Furthermore, the optimal controls are characterized as

The adjoint functions defined in the system (4.6) are bounded.

Proof. The adjoint system in equation (4.6) is linear in the adjoint variables (). Since it is a linear system in finite time with bounded coefficients, it follows that adjoint variables () are uniformly-bounded. The system consisting of the state system, adjoint system and optimal control characterization (equations (4.2), (4.6), (4.8)) is called the optimality system. Since the state and adjoint functions are bounded (from Theorems 4.1 and 4.4), the following theorem that characterizes uniqueness of the optimality system for small time is established. This type of “small time” uniqueness result is standard in nonlinear systems with opposite orientation (Fister, Lenhart, & McNally, 1998).

For a short time period, the solution to the optimality system (equations (4.2), (4.6), (4.8)) is unique.

Numerical simulations of optimal control problem

The optimality system (consisting of equations (4.2), (4.6), (4.8)) is solved numerically using the forward-backward sweep numerical method (Lenhart & Workman, 2007; Mohammed-Awel et al., 2017; Mohammed-Awel & Numfor, 2017). The method requires an initial value on the optimal control pair . The initial condition in (4.3) is used to solve the state system forward in time using the MATLAB built-in ODE45 routine. Next, the adjoint system is solved using the transversality condition (4.7) and the approximated solution of the state system. Then, the value control variables are computed using the control characterization in (4.8), and the control pairs are updated via a convex combination of previous and current values of the control characterization. The process continues until the state variables, adjoint and control values converge.

The parameter values in Table 2 (relevant to the demography and malaria transmission dynamics in Ethiopia (CIA, 2017; Deribew et al., 2017)), and the estimated values of weight constants in the objective functional, are used in the numerical computations. Since the recent estimate for the cost of an ITN is between and (UNICEF), we set . Further, based on the estimations in (White, Conteh, Cibulskis & Ghani, 2011), we set and . Furthermore, and are set at coverage based on the target coverage of IRS and ITNs used in (Griffin et al., 2010). Moreover, the host populations and control functions in the objective functional (4.4) are balanced by choosing appropriate constant weights , , and . Furthermore, it is assumed that in malaria-endemic communities, the order of priorities are (i) to minimize disease in the human population, (ii) to minimize disease in the mosquito population, and (iii) to minimize insecticide resistance in the mosquito population. Thus, the weights , and are chosen based on the ordering . A complete list of the estimated values of the optimal control-related parameters is given in Table 5.

Table 5

Values of the weight coefficients in the objective functional (4.4).

Parameter	Value	Parameter	Value	Parameter	Value	Parameter	Value
A_1	1×10^3	B_1	1.65	C_1	8×10^8	b^max	0.8
A_2	A_1/5	B_2	2B_1	C_2	2C_1	uimax	0.8
A_3	3A_2/4

The population-level effectiveness of the following control strategies will be evaluated:

Strategy 1: Optimal ITNs-only strategy (i.e., and ).

Strategy 2: Optimal IRS-only strategy (i.e., and ).

Strategy 3: Combined optimal ITNs-IRS strategy (i.e., and ).

For each of the aforementioned strategies, the optimal solution for the model (4.2) is solved over a time period (T) of fifteen years (i.e., years) for malaria transmission setting in Ethiopia (unknown parameters are estimated by fitting model (2.4) to data from Ethiopia). Based on the report in (Deribew et al., 2017) (for malaria cases and the estimated Ethiopian population for the year 2000), the following initial values are chosen: , , , , , , , , , , , , , , and . Parameter values for these simulations are as given in Table 2, Table 5, and in the absence of control ITNs and IRS coverages are set at and .

It is, first of all, assumed that chemical insecticides have been used in Ethiopia in the past, and that insecticides resistant vectors exist in the community (so that non-zero initial values of the resistant-type mosquito population can be used in the simulations). Plots for the proportion of infectious humans (that is, the ratio of infectious humans at time t, that is , in the community) are generated. Furthermore, for comparison purposes, the plots for the worst-case scenario (i.e., in the absence of any intervention/control) are generated (and used in these plots). The frequency of resistance allele (denoted by p) is defined as the proportion of resistant-type mosquitoes (that is, ), and the frequency of wild-type allele (denoted by q) is defined as the proportion of wild-type mosquitoes in the mosquito population (that is, ).

Simulations for strategy 1: optimal ITNs-only strategy (that is, and )

The model (2.4) is now simulated to assess the community-wide impact of the singular implementation of an ITNs-only strategy (Strategy 1), and the results obtained are depicted in Fig. 3. It follows from this figure that, in the presence of the optimal level of the ITNs-only control, the frequency of the resistant allele remained at zero for the first 1.5 years, and then increased slowly to 0.5 during the next year. It then increased rapidly to unity a year later (and remained at unity for the remaining duration of the control period) (Fig. 3(a)). Similarly, the frequency of the wild-type allele remained at unity for the first 1.5 years, and then decreased slowly to 0.5 a year later. It then decreased rapidly to zero during the next one year and remained at zero for the remaining duration of the control period (Fig. 3(a)). This figure further shows that, for the worst-case scenario with no intervention (i.e., ), the percentage of infectious humans increased from the initial 8.6%–53% for the first one year, and then decreased to 47% during the next 3 years (and remained there for the next 11 years) (Fig. 3 (b)). On the other hand, in the presence of the optimal level of the ITNs-only strategy, the percentage of infectious humans decreased from the initial to during the first 3 years, and then increased slowly to 4.1% during the next 9 years. It further decreased slowly to 3.3% during the next 3 years (Fig. 3 (b))). Simulations for the profile of infectious humans, in the presence and absence of this ITNs-only control, are depicted in Fig. 3 (c), from which it follows that the implementation of the optimal ITNs-only strategy dramatically decreases the disease burden, in comparison to the worst-case scenario. The profile of the control , depicted in Fig. 3 (d), shows that the optimal ITNs-only control strategy is at the maximum coverage () (for the first 1.2 years), and decreased to 72% during the next 1.6 years. It then increased to 80% during the next 2.3 years and remained at 80% during the rest of the control implementation period. The implication of this result is that the use of the ITNs-only strategy (with coverage at 80%) can significantly decrease the number of infectious humans (from 8.6% to as low as 3.1%) within 3 years. In summary, the community-wide implementation of an optimal ITNs-only strategy (with coverage at 80%) can significantly reduce disease burden (as measured in terms of reduction of the percentage of infectious humans in the population) and effectively manage insecticide resistance during the first 3.5 years of implementation. Unfortunately, however, such insecticide resistance is not effectively managed after the first four years of the 15-year implementation period.

Fig. 3

Numerical simulations of the model (4.2) for Strategy 1. Parameter values used (other than the control parameters b and ) are as given in Table 2, Table 5.

Simulations for strategy 2: optimal IRS-only strategy (that is, and )

The model (2.4) is now simulated to assess the impact of the singular implementation of IRS-only strategy (Strategy 2). The simulations results obtained show that the dynamics of the frequency of both the resistance (p) and wild (q) alleles are similar to the corresponding dynamics observed when Strategy 1 was implemented (see Fig. 3, Fig. 4). Further, when the optimal IRS-only Strategy is implemented, it was seen that the percentage of infectious humans increased rapidly (from the initial to ) for the first 4 years, and then decreased slowly to 40% during the next 6 years (and remained there for the remaining duration of the control period) (Fig. 4 (b)). Furthermore, the results depicted in Fig. 4 (c) show that the disease burden was not reduced with the implementation of an optimal IRS-only strategy. The control profile (depicted in Fig. 4 (d)) shows that the optimal IRS coverage increased quickly to 60% during the first two months, and continued to increase to 80% during the next 10 months (and remained there for the rest of the implementation of the control period). In summary, the community-wide implementation of the optimal IRS-only strategy resulted in an increase in the percentage of infectious humans (from the initial 8.6% to as high as 43% during the first four years of the implementation period) and insecticide resistance (with all adult female vectors becoming wholly resistant after 4 years of implementation). Thus, unlike in the case where an optimal ITNs-only strategy is implemented, the singular implementation of an optimal IRS-only strategy does not offer any population-level benefit, vis a vis decreasing disease burden or effectively managing vector insecticide resistance (in fact, it is detrimental in the long run).

Fig. 4

Numerical simulations of the model (4.2) for Strategy 2. Parameter values (other than the control parameters b and ) are as given in Table 2, Table 5.

Simulations for strategy 3: combined optimal ITNs-IRS strategy (that is, and )

Strategy 3 entails the implementation of the combined optimal ITNs and IRS strategies. Here, too, the dynamics of the allele frequencies p and q are similar to those observed when Strategy 1 or Strategy 2 was singularly implemented (see Fig. 3, Fig. 4, Fig. 5(a)). This figure further shows that this strategy manages insecticide resistance effectively during the first 3.5 years of implementation of the control, and fails to do so thereafter. Fig. 5 (b) shows that the percentage of infectious humans decreased from the initial to during the first 2.5 years, and then decreased slowly to less than 1% during the rest of the control implementation period (Fig. 5 (b)). Further, Fig. 5 (c) shows a dramatic reduction in the number of infectious humans (Fig. 5(b) and (c) show effective disease control, or elimination, after about 10 years of the implementation of this strategy). The profiles of the controls , depicted in Fig. 5 (d) and (e), show that the optimal ITNs coverage is maximum at (for the first 1 year), and decreased to 72.8% during the next 1.6 years, and then increased to (and settled at) 80% during the next 2.3 years. On the other hand, the optimal IRS coverage is initially zero for the first 2 years, and then increased to 50% during the next 3 years. It then decreased slowly to 48% for the remaining duration of the control implementation period. In summary, while the community-wide implementation of the optimal combined ITNs-IRS strategy (Strategy 3), with high enough ITNs coverage, can lead to the elimination of the disease, it fails to effectively manage insecticide resistance (with the frequency of the resistant allele reaching 100% after 4 years of implementation of this strategy).

Fig. 5

Numerical simulations of the model (4.2) for Strategy 3. Parameter values (other than the control parameters b and ) are as given in Table 2, Table 5.

Simulations for strategy 3 for decreased birth rate of new adult resistant-type mosquitoes

In order to assess the impact of the expected fecundity-related fitness cost of insecticide resistance (Alout et al., 2016; Alout et al., 2017), the model (2.4), subject to Strategy 3, is simulated for the special case where the birth rate of new adult resistant-type female mosquitoes () is reduced from the baseline (fitted) value of to (that is, the birth rate is reduced by to account for fitness cost of resistance with respect to fecundity). The simulation results generated for this special case of Strategy 3 are depicted in Fig. 6, from which it follows that the frequency of resistant allele remained at zero for the first 6 years, and then increased slowly to 0.5 during the next 2 years. It then increased rapidly to unity during the next 2 years and remained at unity for the remaining duration of the control period (Fig. 6 (a)). Similarly, the frequency of wild-type allele remained at unity initially (for the first 6 years), and then decreased slowly to 0.5 during the next 2 years. It further decreased to zero during the next 2 years and remained at zero for the remaining duration of the control period (Fig. 6 (a)). Further, the population of infectious humans decreased from the initial to during the first 5 years, and then increased to 3.6% during the next 2 years. It then decreased slowly to less than 1% for the next 3 years (and remained at less than 1% for the remaining duration of the control period) (Fig. 6 (b)). Fig. 6 (c) shows a dramatic reduction in the number of infectious humans, in comparison to the worst-case scenario. Additionally, the profiles of the controls , depicted in Fig. 6 (d) and (e), show that the optimal ITNs coverage rapidly increased from zero to 80% and remained at 80% (for the first one year), then decreased to 46% during the next 4 years. It then increased to 80% during the next 3 years and remained there for the rest of the control period (Fig. 6 (d)). Similarly, the optimal IRS coverage is initially zero for the first 4.5 years, and then increased slowly to 41% during the next 5.5 years and remained there for the rest of the control period (Fig. 6 (e)). In summary, these simulations show that, for the special case of the model (2.4) with the birth rate of resistant mosquitoes is marginally reduced (from to ) to account for the fitness cost of resistance (with regards to fecundity), the community-wide implementation of the combined ITNs-IRS strategy (i.e., Strategy 3) can lead to the effective control of the disease, while also effectively managing insecticide resistance, during the first 8 years of implementation. Unfortunately, however, insecticide resistance develops during the next 7 years of the control implementation period.

Fig. 6

Numerical simulations of the model (4.2) for Strategy 3 where is decreases to from the fitted value , other parameter values (other than the control parameters b and ) are as given in Table 2, Table 5.

Discussion and conclusions

A new mathematical model for malaria transmission dynamics, which incorporates the dynamics of resistant-type and wild-type mosquitoes to insecticides, is presented. Indoor and outdoor mosquitoes are stratified in separate classes, and ITNs and IRS are applied indoors. Back-and-forth mobility of mosquitoes, between indoors and outdoors, is incorporated. The model is rigorously analysed to gain insight into its qualitative features. It is shown that its nontrivial disease-free equilibrium is locally-asymptotically stable when ever a certain epidemiological threshold (denoted by ) is less than unity. The epidemiological consequence of this (local asymptotic stability) result is that, for the case with small number of infectives (in the basin of the disease-free equilibrium) the use of the two control measures (ITNs and IRS) could lead to the effective control (or elimination) of the disease in the community if the epidemiological threshold () can be brought to (and maintain at) a value of less than unity.

Furthermore, optimal control analysis is used to assess the population-level impact of the two vector control strategies considered in this study (i.e., the singular and combined use of IRS and ITNs). In particular, an optimal control problem, based on minimizing infection in the hosts (humans and mosquitoes), insecticide resistance in the vector and the cost of implementing the two controls in the community, is formulated (and solved numerically). Numerical simulations of the resulting optimal control problem, using parameter values obtained from fitting the model (2.4) with data relevant to malaria transmission dynamics in Ethiopia, show that:

the singular implementation of an optimal ITNs-only strategy (Strategy 1) in the community could lead to a dramatic reduction in the disease burden in humans (e.g., it reduces the percentage of infectious humans from the initial assume 8.6% to as low as 3.3%) and effectively manage insecticide resistance during the first 3.5 years of its implementation. Unfortunately, however, it fails to continue to effectively manage the insecticide resistance after the first 4 years (Fig. 3);

unlike for the case for the singular implementation of an optimal ITNs-only strategy, the singular implementation of an optimal IRS-only strategy (Strategy 2) in the community fails to lead to either the effective control of the disease or manage insecticide resistance. In fact, it resulted in an increase in the percentage of infectious humans in the community from the initial 8.6% to as high as 40%. Furthermore, it fails to manage insecticide resistance after the first 4 years of its implementation (Fig. 4);

like for the case of the singular implementation of an optimal ITNs-only strategy, the implementation of the combined optimal ITNs-IRS strategy (Strategy 3) in the community resulted (expectedly) in a dramatic reduction in the disease burden (in particular, it reduces the percentage of infectious humans in the community from the initial 8.6% to as low as 1%) (Fig. 5), in addition to effectively managing insecticide resistance during the first 3.5 years of its implementation (albeit is fails to manage such insecticide resistance thereafter);

for the case when the expected fitness cost of insecticide resistance with respect to fecundity is accounted for (by reducing the birth rate of new adult resistant-type female mosquitoes () from the fitted value of to ), the combined optimal ITNs-IRS strategy (Strategy 3) led to the effective control of the disease, as well as manage resistance effectively during the first 8 years of implementation of the strategy (this strategy fails to manage such insecticide resistance thereafter) (Fig. 6).

This study shows that the singular use of an optimal ITNs-only strategy, or in combination with optimal IRS, can lead to effective control of the disease, while also effectively managing insecticide resistance during the first few years of its implementation, in the malaria-endemic setting (Ethiopia) considered. It is worth emphasizing that one possible explanation for the fact that some of our results (Fig. 3, Fig. 5) show effective disease control but insecticide resistance is not effectively managed, is that the fitness costs associated with insecticide resistance (measured by the fitted parameter values) is high. When these costs are reduced (e.g., when the birth rate of resistant mosquitoes is reduced, to account for the fitness cost of resistance with respect t fecundity), the optimal combined ITNs-IRS strategy controlled both the disease and insecticide resistance for the first 8 years (Fig. 6). This result is in agreement with the results reported in (Barbosa & Hastings, 2012; Brown et al., 2013) that fitness costs are key elements in disease control and insecticide resistance management strategies. Our study, therefore, calls for further lab experimentation by entomologists (and modeling work) to obtain data that can be used to obtain improved estimates of the fitness-related parameters. This is crucially needed for gaining better understanding (and realistically quantifying) the population-level impact of insecticide resistance on the epidemiology and control of malaria. In summary, this study shows that the prospect of the effective control of malaria spread (while minimizing the risk of insecticide resistance in the female adult mosquito population), using ITNs alone or its combination with and IRS, are promising, provided that the effectiveness and coverage levels of these (ITNs-only or the combined ITNs-IRS) interventions are at optimal levels.