Skip to main content

Optimal control of malaria: combining vector interventions and drug therapies

Abstract

Background

The sterile insect technique and transgenic equivalents are considered promising tools for controlling vector-borne disease in an age of increasing insecticide and drug-resistance. Combining vector interventions with artemisinin-based therapies may achieve the twin goals of suppressing malaria endemicity while managing artemisinin resistance. While the cost-effectiveness of these controls has been investigated independently, their combined usage has not been dynamically optimized in response to ecological and epidemiological processes.

Results

An optimal control framework based on coupled models of mosquito population dynamics and malaria epidemiology is used to investigate the cost-effectiveness of combining vector control with drug therapies in homogeneous environments with and without vector migration. The costs of endemic malaria are weighed against the costs of administering artemisinin therapies and releasing modified mosquitoes using various cost structures. Larval density dependence is shown to reduce the cost-effectiveness of conventional sterile insect releases compared with transgenic mosquitoes with a late-acting lethal gene. Using drug treatments can reduce the critical vector control release ratio necessary to cause disease fadeout.

Conclusions

Combining vector control and drug therapies is the most effective and efficient use of resources, and using optimized implementation strategies can substantially reduce costs.

Background

Vector-borne diseases inflict significant levels of human morbidity and mortality. Current estimates suggest that vector borne disease account for 17% of the global disease burden and over half of the world’s population are at risk of contracting a vector-borne disease. Just under half the world’s population live at risk of dengue [1] with more than 300M people contracting dengue annually [1]. Similarly, while the global incidence of malaria is falling, in 2015,  200M malaria infections were reported leading to an estimated 429,000 deaths mostly among African children under the age of 5 [2].

The World Health Organization (WHO) recommends insecticide-treated bed nets (ITNs), indoor residual spraying (IRS) and anti-malarial drug therapies, specifically, the treatment of clinical malaria with artemisinin-based combination therapy (ACT), as the principal methods used to combat malaria [2, 3]. Catalyzed by the Roll Back Malaria Initiative around the United National Millennium Development Goals, a widespread scale-up of coverage of these control interventions successfully reduced and locally eliminated malaria in sub-Saharan Africa; between 2000 and 2015, Plasmodium infection in endemic regions of Africa halved and the incidence of clinical disease fell by 40% [4]. This remarkable and widespread reduction is estimated to have averted 663 million clinical cases of malaria since 2000 [4]. However, these gains are fragile and they are increasingly threatened by the emergence of Plasmodium strains that are resistant to anti-malarial drugs and mosquitoes that are resistant to the insecticides used to kill them.

Due to widespread resistance to anti-malarial drugs, the WHO recommends ACT as the front-line treatment in all countries with endemic malaria. ACT was first introduced in the mid-1990s and are generally cost-effective [5]. Artemisinin is a highly potent plant-based compound which clears parasitaemia more rapidly than all other currently available anti-malarial drugs [6, 7]. Moreover, genetic mutations which are thought to confer resistance to the anti-malarial drug chloroquine also increase susceptibility to artemisinin and quinine [8]. Therefore, ACT can help delay the development of resistance, as they combine fast-acting artemisinin with another class of anti-malarial drug (such as quinines or anti-folates). Despite these efforts, artemisinin resistance first emerged in the mid-2000s in Cambodia, resulting in longer parasite-clearance times and rising failure rates for some artemisinin-based drug combinations [9]. This resistance has now emerged in or spread across mainland Southeast Asia, appearing in Thailand, Cambodia, Myanmar and Vietnam [10,11,12]. For now, ACT is still an effective treatment in these countries, as parasites remain susceptible to some partner drugs [13]. However, it is only a matter of time before these drugs begin to fail and while new anti-malarial drugs are being developed, it will be several years before they become available.

By far the most cost-effective intervention for reducing Plasmodium infection rates since 2000 has been the widespread roll-out of ITNs—bed nets treated with pyrethroid insecticides [4, 14]. As ITNs kill or disable the mosquitoes which land on them, even modest adoption rates of ITNs can reduce the malaria vector population and achieve community-wide benefits [15]. Unfortunately, although not surprisingly, the intense exposure to insecticides due to the adoption of ITNs and increase in IRS programmes is driving the spread of insecticide resistance in mosquitoes. Pyrethroid resistance was first detected in the mid-1990s [16], and is now ubiquitous in all African malaria vectors and is increasing in strength, meaning mosquitoes can tolerate ever-higher levels of chemical exposure [17]. Several countries have identified mosquito populations that are starting to develop resistance to all four classes of insecticide that are approved by the WHO for IRS (organochlorines, organophosphates, carbamates and pyrethroids [18,19,20]). Worryingly, there are indications that this insecticide resistance is already compromising the effectiveness of ITN and IRS control measures [17].

In this environment of ever-growing resistance, novel vector control techniques which seek to reduce mosquito populations without the need for chemical insecticides may be vital for effective malaria management [10]. Pest control via the release of sterile insects is an old idea [21] and its successes have been well-documented in other species [22, 23]—although unforeseen circumstances, such as unidentified wild breeding sites, can cause problems for these control programmes [24, 25]. The Sterile Insect Technique (SIT) may offer a means to control mosquito populations [26], and much modelling work has shown that SIT and its transgenic equivalents could be successful [27,28,29,30]. The traditional radiation-based SIT represents a form of ‘early-acting lethality’, in the sense that offspring die at the very earliest stage. However, such early-acting lethality may be undesirable: ecological studies show that larval habitats can exhibit overcompensatory density dependence, where the size of the adult population increases as larval density decreases [31]. Consequently, traditional SIT control may be ineffective [32] or, worse, even increase adult mosquito populations [33, 34]! If overcompensatory density dependence occurs, an SIT is needed which exhibits ‘late-acting lethality’, where larvae fully contribute to density-dependent competition but then die at the end of the larval stage. Advances in molecular genetics are making such late-acting SITs possible as transgenic insects carrying a dominant lethal gene have been developed [32, 35, 36], where the progeny of transgenic mosquitoes die at the end of the larval or pupal stage [37].

Information on the costs and effectiveness of different interventions are essential so that policy makers can make informed decisions. Detailed studies of the efficacy and cost-effectiveness of malaria interventions show that ACT is generally cost-effective [5, 38], and ITNs and IRSs are highly cost-effective [39,40,41]—though in the long term resistance is likely to render them ineffective. Therefore, even if it is technically possible to engineer transgenic insects which could reduce vector populations to a level where they can no longer sustain malaria transmission, if they are prohibitively expensive compared to the alternatives they are unlikely to change malaria management policy. For this reason, it is important to assess the cost-effectiveness of novel vector control techniques at the earliest opportunity; although such estimates require many assumptions to be made, they can still be indicative of whether policy makers should take such vector control techniques seriously or not. To date, little is known about the the cost-effectiveness of early- or late-acting SITs. Cost–benefit analyses have been performed for traditional SITs in other pest insect species [42], but not for mosquitoes. One study has examined the cost-effectiveness of using constant releases (rather than dynamically managed releases, which is investigated here) of transgenic late-acting SIT mosquitoes to reduce dengue [43]. However the cost-effectiveness of dynamic, optimized releases of transgenic late-acting SITs for malaria has not been examined, nor has any comparison been made between the cost effectiveness of early vs late-acting SIT.

Here, conventional early-acting SIT with genetics-based late-acting SIT is compared for various strengths of larval density dependence to assess how ecology affects the efficacy and cost-effectiveness of vector control. This vector control model is coupled to an epidemiological model of malaria spread between vector and host, using P. falciparum estimates for parameterization, under the action of ACT. The aim is to investigate both the efficacy of the combined drug and vector control strategies and identify the most cost-efficient way to deploy these strategies. To do this, two control parameters are introduced: the SIT release ratio u and artemisinin treatment ratio w, which act on the population model and disease model, respectively. Cost functions for these two parameters are developed and then weigh against the economic burden of endemic disease. This ‘optimal control framework’ is used to investigate the most cost effective way to reduce the disease burden.

Methods

A stage-structured model of the mosquito population [44, 45] is combined with a Ross–MacDonald epidemiological model of the malaria dynamics [46,47,48,49,50]. The governing ordinary differential equations for the larval (L) and adult female (A) mosquito populations and the proportion of the human (h) and vector (v) populations that are infected (and infective) are

$$\begin{aligned} \frac{\mathrm {d}L}{\mathrm {d}t} =&\, \rho A(t) C(A,u) - f(L) L(t) - (m + \mu _L) L(t) ,\end{aligned}$$
(1a)
$$\begin{aligned} \frac{\mathrm {d}A}{\mathrm {d}t} =&\, \frac{m}{2} L(t) - \mu _A A(t), \end{aligned}$$
(1b)
$$\begin{aligned} \frac{\mathrm {d}h}{\mathrm {d}t} =&\, \frac{A(t)}{N} b a (1 - h(t)) v(t) - (\gamma + w(t) s) h(t), \end{aligned}$$
(1c)
$$\begin{aligned} \frac{\mathrm {d}v}{\mathrm {d}t} =&\, b c (1 - v(t)) h(t) - \left( \mu _A + \frac{1}{A(t)}\frac{\mathrm {d}A}{\mathrm {d}t}\right) v(t), \end{aligned}$$
(1d)

where \(\rho\) is the mosquito oviposition rate, m is the rate at which larvae mature into adults (m/2 represents half the larvae maturing into adult females), \(\mu _L\) and \(\mu _A\) are the density-independent mortality rates of the larval and adult stages, respectively, N is the total host population, b is the biting rate of the vectors, a is the proportion of bites by infected mosquitoes on susceptible humans that produce an infected human, c is the proportion of bites on infected humans by a susceptible mosquitoes that produce an infected mosquito, \(\gamma\) is the rate of recovery. The disease model is parameterized using P. falciparum estimates, see Table 1. ACT is assumed to affect a fraction w (the treatment proportion, or drug coverage) of infected humans through an increase in their recovery rate by an amount s. A continuous-time approach is chosen over a discrete-time model in order to capture the effect of control technologies and density-dependent regulation in overlapping generations. However, in order to develop the mathematics for an optimization analysis, a simplifying assumption is made such that no time-delayed effects occur, such as parasite incubation, egg development or drug-clearing times.

The function C(Au) in (1) describes the effective reduction in fecundity due to mating with released sterile males (assuming the populations are well-mixed and undergo random mating), and is defined as

$$\begin{aligned} C(A,u) = \frac{A(t)}{A(t) + u(t)A^{*}}, \end{aligned}$$
(2)

where \(A^{*}\) is the equilibrium adult female mosquito population (without vector control), and u is the release ratio of sterile male mosquitoes (hence \(u A^{*}\) gives the instantaneous, i.e. daily, number of insects being released). The equilibrium levels of endemic P. falciparum infection (as a fraction of the total populations) are \(h^*\) (hosts) and \(v^*\) (vectors). Using Global Health Observatory data [51] for malaria prevalence in the most hard-hit African nations and estimates of sporozoite rates from Kilama et al. [52] to produce order-of-magnitude approximations for \(h^*\) and \(v^*\), we use the transmission efficiencies a and c to set the endemic equilibrium levels to \(\sim\)1%.

While there is some rudimentary information on the magnitude and type of intraspecific competition in larval Anopheles [31, 53], this is neither sufficiently comprehensive nor definitive. As such a relatively flexible form of density dependence [54] is chosen that has been recently used elsewhere to explore the impact of ecological intraspecific competition and genetics-based methods of vector control [34]. This form of density dependence was originally proposed by Maynard Smith and Slatkin [55] and is of the form

$$\begin{aligned} f(L) = \ln {\left( 1 + (\nu L)^{\beta }\right) }. \end{aligned}$$
(3)

Density dependent effects are set by a scale parameter (\(\nu\)) and all larval cohorts experience the same strength of density dependence (which can range from contest—resources monopolized by few individuals—to scramble—resources shared equally by all—by varying the coefficient \(\beta\)). Both the timing of density dependence with respect to genetic lethality mechanisms and the strength of density dependence are well-known to affect vector control programmes and can lead to ‘bad’ SIT effects [33, 34]. \(\nu\) is set to satisfy the relation

$$\begin{aligned} \nu = \frac{m}{2 k^* N \mu _A} \left( \mathrm {e}^{\{\rho m/(2 \mu _A) - m - \mu _L\}} - 1 \right) ^{\frac{1}{\beta }}, \end{aligned}$$
(4)

to ensure that the equilibrium population \(A^*\) is equivalent to \(k^* N\), the number of vectors per host at equilibrium multiplied by the host population.

The population model (1a) and (1b) is modified slightly if sterile insect releases are replaced with releases of insects modified with a late-acting lethal gene [an early-acting lethal gene is adequately modelled by (1a) and (1b)]. Genetic mortality is assumed to occur at the end of the larval stage, such that all offspring will contribute fully to density-dependent mortality. Thus, the stage-structured mosquito population model takes the form

$$\begin{aligned} \frac{\mathrm {d}L}{\mathrm {d}t} =&\, \rho A(t) - f(L) L(t) - (m + \mu _L) L(t) \end{aligned}$$
(5a)
$$\begin{aligned} \frac{\mathrm {d}A}{\mathrm {d}t} =&\, \frac{m}{2} L(t) C(A,u) - \mu _A A(t), \end{aligned}$$
(5b)

where C(Au) is as defined in (2). Model (5) will be referred to as ‘late-acting SIT’ as opposed to the ‘early-acting SIT’ (Eq. 1).

Table 1 Parameter definitions and values. Disease parameters use P. falciparum estimates

Optimal control

A cost functional is defined that measures the cost per capita, in US$, of disease level (h), artemisinin medication (w) and vector control (u) over a fixed time period \(t\in [0,T]\):

$$\begin{aligned} J_{pc}[\varvec{x}(u,w)] = \int _0^T \mathrm {e}^{-\psi t} \left[ \theta _{1} h^{m_{1}} + \theta _{2} h w^{m_{2}} + \theta _{3} k^{*} u^{m_{3}}\right] \mathrm {d}t, \end{aligned}$$
(6)

where future costs are discounted at a constant rate \(\psi\) to deal with opportunity costs in investing in public health control for immediate benefits [56, 57], \(\varvec{x} = (L, A, h, v)^{\top }\) is the vector of state variables and the parameters \(\theta _i\), \(k^*\) and \(\psi\) take values given in Table 1. The price factors \(\theta _i\) and the cost function exponents \(m_i\) govern the scale and behaviour of cost for each contributing factor, and form the (per capita) cost functions \(C_h = \theta _{1} h^{m_{1}}\), \(C_w = \theta _{2} h w^{m_{2}}\) and \(C_u = \theta _{3} k^{*} u^{m_{3}}\) (see Appendix A: “Optimal control formulation” for more details). Drug coverage was capped at a maximum of \(60\%\) (\(w_{\mathrm {max}} = 0.6\)); this prevents the optimization approach from recommending drug coverages close to \(100\%\), as such high levels would be unachievable in reality.

The goal is to minimize the total cost of the malaria management strategy, where the controls are the insect release ratio (u) and the artemisinin treatment ratio (i.e. drug coverage, w):

$$\begin{aligned} \mathop {\mathop {\min }\limits _{u \ge 0} }\limits _{0 \le w \le 0.6} \{ {J_{pc}}\}. \end{aligned}$$
(7)

The method used to solve this optimisation problem is described in Appendix A: “Optimal control formulation”, along with some associated mathematical results.

Results

Disease elimination and extinction

Model simulations (based on the parameters in Table 1) show, as expected, that in the absence of control, the proportion of infected individuals, h, and infected vectors, v, will reach an endemic level (Fig. 1a). Introducing drug control, whereby individuals recover at an enhanced rate, can reduce endemic levels of disease in hosts (Fig. 1b). Mosquito vector control (through population suppression technologies such as SIT or genetic engineering) can effectively reduce disease (Fig. 1c); by combining both vector control and drug-based therapies, disease spread is reduced and disease control is most effective (Fig. 1d).

Fig. 1
figure 1

Vector control and drug therapies are most effective when used in tandem. Proportions of humans (h) and vectors (v) infected are plotted under four control scenarios. a no control, \(w=0\), \(u=0\); b only artemisinin treatment, \(w=0.05\) (\(5\%\) drug coverage), \(u=0\); c only vector control, \(w=0\), \(u=0.2\) (releasing modified males at a rate of \(20\%\) of the wild male populations per day); d both artemisinin treatment and vector control, \(w=0.05\), \(u=0.2\). The total number of infected mosquitoes at time \(t=90\) days for each scenario is a 6125 b 1145 c 238 d 55; the vector control suppresses the vector population significantly. Early-acting SIT is assumed

Fig. 2
figure 2

Vector control broadens disease-free parameter space, causing disease fade-out even for high biting rates. The level of endemic disease, \(h^*\), across bw parameter space (biting rate–drug treatment proportion) with a no vector control, b early-acting SIT releases and c late-acting SIT releases. The release ratio in b and c is \(u=0.1\) (daily releases of \(10\%\) of the wild male population), with the ratio being of the current vector population A(t) rather than the equilibrium population \(A^*\) (elsewhere in the paper the ratio of \(A^*\) is used). Disease-free regions of parameter space are coloured white. Other parameters as in Table 1

Fig. 3
figure 3

Vector control releases can achieve disease fade-out for sparse drug coverage when biting rates are below \(b<1.2\). The critical release ratio \(u=u_{\mathrm {c}}\) that leads to disease fadeout (\(h^*\rightarrow 0\)) for a given treatment proportion w as the mosquito biting rate b varies, as calculated using (46) and (47). The release ratio u is of the current vector population A(t) rather than the equilibrium population \(A^*\) (elsewhere in the paper the ratio of \(A^*\) is used). Early and late SIT are compared; all other parameters as in Table 1

It is usual to discuss values of the basic reproductive ratio (\(R_0\)) of disease when analysing epidemiological models. The value of \(R_0\) is a static measure of the speed of disease spread in a naïve environment [49, 58]. Effective reproductive ratios can be constructed which take account of growing host populations [59] or heterogenous environments [60], but these still cannot capture the dynamic effect of vector control via SIT or transgenic insect releases. Rather than investigating how quickly a disease might invade a disease-free population (a situation in which SIT releases would not be in use), here it is investigated how far below the control-free endemic equilibrium the prevalence of disease is pushed by vector control (a situation in which insect releases would be called for). The equilibria of the system with and without control are given in Appendix B: “\(R_0\) and equilibria”. If the application of control (drug-based and vector) takes the new equilibrium of h to zero, then that point in parameter space is considered disease-free. In the classic control-free, disease-naïve scenario, the \(R_0\) value depends upon the square of biting rate b, and is linear in a and c (see Appendix B: “\(R_0\) and equilibria”). Thus, b is of greater importance to disease spread than the other parameters that constitute \(R_0\) (a, c, k, \(\gamma\) and \(\mu _A\)). Changes in b will therefore exert relatively large changes in disease prevalence, which can be balanced by implementation of drug therapies (w) and vector control. Therefore, disease-free regions that occur in bw parameter space are investigated.

For the constant release ratio \(u=0.25\) (where here releases are made in a ratio to the current vector population A(t) rather than the equilibrium population \(A^*\), a comparison first made in [37]) the disease-free region of bw space is increased considerably (Fig. 2), and the severity of the endemic disease away from the disease-free region is reduced for both early- and late-acting SIT (compare Fig. 2a with Fig. 2b, c). The greatest effect of vector control on \(h^*\) is seen at \(w=1\) (Fig. 2b, c), showing that the most effective way to combat endemic disease is through a combination of vector control and drug-based therapies. Late-acting SIT is seen to out-perform early-acting SIT at the same release ratio.

The critical release ratio \(u = u_{\mathrm {c}}\) (when releasing uA(t) insects rather than \(u A^*\) insects) that leads to disease fadeout (\(h^*\rightarrow 0\)) may be found (see Appendix B.2: “Equilibria and the effect of control”). For biting rates below \(b\approx 1\) (i.e. vectors blood feed as frequently as once per day), manageable vector control releases, \(u<0.1\), can be used to cause disease fadeout when drug coverage is sparse, \(w<0.2\) (Fig. 3). For high biting rates, \(b>1.2\), the use of drug therapies can substantially reduce the critical release ratio, for example for \(b=2\) maintaining drug coverage of \(60\%\) can reduce the required mosquito release rates from \(30\%\) of the wild male population daily to under \(5\%\).

Optimal control

Control regimes may more efficiently (with respect to time and cost) combat disease when the optimal control solutions (described in Appendix A.1: “Optimality conditions”) are implemented (Fig. 4). In Fig. 4a the four naïve control scenarios of Fig. 1 are compared with an optimal vector control strategy \(u^*(t)\), (29), with no artemisinin treatment, and an optimal control strategy using both vector control and drug therapies \(w^*(t)\), (18). The breakdown of costs (using the pricing parameters given in Table 1) neglects the one-off capital investment used to a construct insect-rearing facility (the facility could be used again for future control efforts—the operational costs for a single season are presented here).

When density-dependent effects at the larval stage are weak, early- and late-acting SIT technologies are similar in efficacy; when density-dependent effects are strong, the two technologies diverge and late-acting SIT is shown to be the most cost-effective (Fig. 5). Although the form of the optimal control strategies for early- and late-acting SIT are similar (Fig. 5b, c), the enhanced efficacy of late-acting SIT due to density-dependent mortality means lower release ratios produce optimal results, requiring facilities with lower weekly rearing capacities.

Fig. 4
figure 4

Optimal strategies for releases and drug treatments can substantially reduce the cost of managing malaria. a The proportion of infected humans under six control scenarios (where w is the host treatment proportion and u is the insect release ratio with respect to the wild vector population): no control (x-axis label “0,0”), \(w=w_0\), \(u=0\) (x-axis label “\(w_0\),0”), \(w=0\), \(u=u_0\) (x-axis label “0,\(u_0\)”), \(w=w_0\), \(u=u_0\) (x-axis label “\(w_0\),\(u_0\)”), \(w=0\), \(u=u^*(t)\) (x-axis label “0,\(u^*\)”), \(w=w^*(t)\), \(u=u^*(t)\) (x-axis label “\(w^*\),\(u^*\)”), where \(w_0 = 0.05\), \(u_0 = 0.2\) (\(5\%\) drug coverage and releasing \(20\%\) of the wild male population per day) and \(w^*\), \(u^*\) are the optimal control strategies defined in (18) and (29), respectively. b The total cost of the scenario, including spending on traditional healthcare (h), spending on artemisinin treatment (w) and spending on insect releases (u). Early-acting SIT is assumed. Quadratic cost functions of h, w and u are assumed

Fig. 5
figure 5

Late-acting lethality suppresses the effects of over-compensatory larval density dependence. a Cost per case averted (excluding initial capital investment in construction) for optimized release strategies of early-acting and late-acting SIT, for the parameters given in Table 1, as the strength of the larval density dependence \(\beta\) is increased. The density dependence runs from contest for \(\beta < 1\) to scramble for \(\beta > 1\). The number of cases averted is approximated by ([mean no. with no control − mean no. with control]/average duration of disease) \(\times\) no. of days). The change in the optimal release strategy is shown for b early SIT, and c late SIT for three specific values of \(\beta\). Artemisinin treatment is not used

Fig. 6
figure 6

Vector control costs, unlike medical care costs, are insensitive to the form of cost function. a Cost per capita for the optimal strategy over a 40 day control period for three cost functions: linear, quadratic and square root; and three epidemiological parameter sets: low, medium and high (for low \(a,b,c,k^* = 75\%\) of medium, for high \(a,b,c,k^* = 125\%\) of medium). In the top figure of a the cost function \(C_h\) is varied while keeping \(C_u\) and \(C_w\) quadratic (\(m_2 = m_3 = 2\)); in the bottom figure the cost function \(C_u\) is varied while keeping \(C_h\) and \(C_w\) quadratic (\(m_1 = m_2 = 2\)). The number of cases averted is approximated by ([mean no. with no control - mean no. with control]/average duration of disease) \(\times\) no. of days). b A sketch of the three cost functions \(\hat{\theta }_i (\cdot )\) (linear), \(\hat{\theta }_i (\cdot )^{1/2}\) (square root) and \(\hat{\theta }_i (\cdot )^2\) (quadratic), where the altered price \(\hat{\theta }_i\) that multiplies the cost function is: 0.8 (square root), 1 (linear), 1.5 (quadratic). Late-acting SIT is used

The general cost function exponents in (6) for h and u allow for possible nonlinearities in cost to be investigated separately for each variable; the exponent of \(C_w\) is restricted to \(m_2>1\) as outlined in the discussion in Appendix A.1: “Optimality conditions”. Figure 6a shows that reducing the cost exponent \(m_1\) of \(C_h\) increases the cost per capita over the control period for the combined vector control and drug therapies optimal strategy. An exponent \(m_1<1\) describes a situation where treating a small proportion of the population is expensive, but the cost levels off at higher treatment proportions as economies of scale take effect; conversely, an exponent \(m_1>1\) describes a situation where treating a reasonable proportion of the population can be done cheaply, but the price steepens for high proportions, capturing the fact that clinics may exceed capacity, or some sections of the population may be more difficult to reach (see Fig. 6b). Changing the cost function of the vector control variable u has a relatively small effect due to the low cost associated with releasing a large number of mosquitoes (Fig. 6a, lower plot). The optimal strategy for vector control is generally to swamp the wild population at the beginning of the control period, and then reduce the releases as the control takes effect (Fig. 7a). For a square root cost function it is very cheap to release a large number of mosquitoes (which is mathematically optimal), but the greater the initial influx of released mosquitoes, the greater the capacity must be of the rearing facility. This can lead to construction costs that dwarf the control-free costs of conventional treatment of endemic malaria over a single season. In fact, diminishing returns are seen in the effect on h of increasing the initial release ratio u (Fig. 7b), as the square root cost function solution for u causes only a slight decrease in h over the control period when compared to linear or quadratic cost functions, for which the optimal release strategies are initially more conservative.

Concentrated or distributed control?

A small network of three equally spaced populations is established to investigate the effects of migration on the efficacy of vector control, Fig. 8 (see Appendix C: “A network model” for model details). Three methods of vector control are compared: a concentrated control approach, where one population is subjected to a constant high SIT release ratio \(u = 3\), a distributed control approach, where each population is subjected to a constant low SIT release ratio \(u=1\) (the total number of released insects is the same in these two cases), and an optimal vector control approach for each population, with and without combined optimized drug treatment. As SIT methods are an area-wide control strategy, mosquito dispersal, although often limited by energy resources, from an uncontrolled population can damage the efficacy of vector control in a focal population (Fig. 8a, b). Distributing control efforts across populations is shown to be a more economically efficient and effective vector intervention strategy than focusing on a single disease hot-spot (Fig. 8b, c). Optimal vector control strategies (rather than constant release ratios) can reduce total expenditure, but the greatest monetary saving is made when combining optimal vector control with optimal drug treatment (Fig. 8d–f).

Discussion and conclusions

Here a novel framework that links disease dynamics, an ecological mosquito population model and an economic optimal control formulation has been developed. An optimal control framework was used to assess the efficacy and cost-effectiveness of combining ACT with two novel vector control technologies: releasing either sterile mosquitoes or transgenic mosquitoes carrying a late-acting lethal gene. In all cases, the transgenic late-acting mosquitoes are more cost-effective than the early-acting sterile insects. This is because the transgenic mosquito larvae compete with the wild mosquito larvae; the resulting competition increases larval mortality which reduces the wild type adult mosquito population. Immigration of vectors from uncontrolled populations can damage a concentrated control effort, and found that it is optimal to use distributed control strategies combined with drug treatments. Results from this analysis revealed that the combined drug and vector control interventions have a cost per case averted of \(\sim \$1\) (2016 US$) if construction costs are neglected, or \(\sim \$4\) if construction costs are included in the control costs for a single season. These estimates are necessarily based on assumptions about the cost of transgenic releases and construction costs of insect rearing facilities, which are uncertain given that this technique is at such an early stage of development. Nevertheless, even if the real costs are double or even triple those predicted here, they compare extremely favourably with the economic costs of the current leading malaria control methods [61].

Boundaries of disease-free regions of (biting rate–drug treatment ratio) parameter space were found by calculating post-control population equilibria. This method was chosen over the usual \(R_0\) calculations for two reasons: (i) the dynamic effect of vector control could not be captured by the static \(R_0\) value, (ii) the disease-naïve invasion basis of \(R_0\) does not fit well with the goal of assessing the control of endemic disease. Vector control has a positive effect on the size of disease-free regions and can reduce the prevalence of disease outside of these regions (Fig. 2). The use of drug therapies can significantly lower the critical release ratio required to eradicate disease even for high vector biting rates (Fig. 3), showing that combining control methods can be an effective disease management strategy.

Drug therapies act to reduce the level of endemicity and suppress the initial spread of disease, while vector control combats disease over time and continuously reduces infection levels towards zero (Fig. 4a). Optimal control strategies follow a robust pattern, indicating that a high initial release ratio is necessary for cost-effective vector control (Fig. 7a). A large proportion of the cost of a single-season control regime is taken up by the construction of the insect rearing facility, but these costs were neglected as it was assumed that vector control efforts will be repeated every season using the same infrastructure. Multi-season control efforts will benefit more from optimal vector release strategies as returns are made from the construction investment in the form of efficiency savings and reductions in the burden of conventional health care. There is an opportunity to extend this work by including construction costs in the optimization procedure which, while not straightforward, would enable an investigation of the (possibly competing) dual desires to keep initial costs down but strive for long-term cost efficiency.

An investigation of nonlinear (concave or convex) cost structures was carried out to model different real-world cost responses (Fig. 6). For concave-cost SIT releases, diminishing returns were found in the effect on the infection level when swamping the wild population with large numbers of lab-reared insects (Fig. 7). This suggests that it is wise to be conservative with maximum facility output requirements and highlights the importance of checking (prior to implementation) the cost-effectiveness of optimal strategies derived using a variety of cost structures. The exact form of the cost functions could be better-informed with additional, or more specific, data. Alternatively, to develop better comparisons between different regions, costs associated with endemic malaria and costs of controls might be weighted using GDP per capita, or through a purchasing power parity scheme (such as the Geary–Khamis dollar; [62]). Vector control costing could be made more reliable if information on Anopheles gambiae facility operational costs (insect rearing, staff pay, heating, materials etc.) was available, rather than inferring from data about facilities rearing different species. Information on the cost of surveying (needed for estimating current vector populations and for effectively targeting releases) would also be beneficial.

Fig. 7
figure 7

Increasing the initial insect release ratio has diminishing returns on the effect on disease prevalence. Example solutions from the default epidemiological parameter set of Fig. 6a when \(C_u\) is varied through square root, linear and quadratic. a Shows the optimal vector release strategy, b shows the resulting suppression of disease in the host population. Late-acting SIT is used. Other parameters as in Table 1

Fig. 8
figure 8

Migration can harm concentrated control efforts; distributed vector control with combined drug treatment is optimal. The proportion of hosts infected in each of three identical populations with SIT release ratios \(u_1\), \(u_2\) and \(u_3\), where migration is \(10\%\) between nodes except in a where there is no migration. a \(u_1=3\), \(u_2 = u_3 = w_i = 0\); b \(u_1=3\), \(u_2 = u_3 = w_i = 0\) (concentrated control); c \(u_1 = u_2 = u_3 = 1\). \(w_i = 0\) (distributed control); d \(u_i(t) = u_i^*(t)\), \(w_i = 0\) (optimal vector control); e \(u_i(t) = u^*(t)\), \(w_i(t) = w_i^*(t)\) (optimal combined drug treatment and vector control). f Cost breakdowns for the five scenarios plotted above. All parameters as in Table 1. Late-acting lethality is used

A network population model with migration between nodes weighted by relative distance allowed optimal control in a spatial environment with identical host populations to be studied. Migration from uncontrolled nodes was seen to lessen the effectiveness of control on neighbouring nodes (Fig. 8a, b). Distributed control at all pertinent nodes was shown to be the most efficient use of resources (Fig. 8c). While it is difficult to parameterize a model of mosquito movement (though values between 4–\(24\%\) migration between neighbouring villages are reported in [63]), the spatial ecology of A. gambiae will undoubtedly be crucial in determining the spread or containment of control efforts. The network model proposed here could be improved by including human movement [64] as a means to spread infection. Using the optimal control framework in such a model would require reformulating cost functions capable of capturing heterogeneous populations, which is beyond the scope of this study.

The ecological and epidemiological models used here are simplified to allow an optimal control framework to be established. Worthwhile extensions to this work would be the addition of a disease incubation period in either the vector or host; use of time delays between egg laying and hatching and pupation and emergence; modelling possible periods of immunity to reinfection conferred by artemisinin treatments (if the prophylactic period is greater than the infectious period; [65]), although the relationship between drug half-life and resistance development is complex [66, 67]; a relaxation of the assumptions of random mating between wild and released mosquitoes (unsuccessful mating has been a serious impediment to codling moth control programmes; [68]) and of perfect sterility and lethality; the inclusion of fitness costs for sterile or transgenic mosquitoes. Finally, that other vectors of malaria (e.g. Anopheles funestus) can be important [69], and interspecific competition may lead to a shift of primary vector if A. gambiae is suppressed [70]. A model incorporating the ecologies of two vector species would shed light on the cost-effectiveness of control if primary vector replacement occurs.

References

  1. Bhatt S, Gething PW, Brady OJ, Messina JP, Farlow AW, Moyes CL, Drake JM, Brownstein JS, Hoen AG, Sankoh O, Myers MF, George DB, Jaenisch T, Wint GRW, Simmons CP, Scott TW, Farrar JJ, Hay SI. The global distribution and burden of dengue. Nature. 2013;496(7446):504–7.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  2. WHO. World malaria report 2016. Geneva: World Health Organization; 2016.

    Google Scholar 

  3. WHO. World malaria report 2013. Geneva: World Health Organization; 2013.

    Google Scholar 

  4. Bhatt S, Weiss DJ, Cameron E, Bisanzio D, Mappin B, Dalrymple U, Battle KE, Moyes CL, Henry A, Eckhoff PA, Wenger EA, Briet O, Penny MA, Smith TA, Bennett A, Yukich J, Eisele TP, Griffin JT, Fergus CA, Lynch M, Lindgren F, Cohen JM, Murray CLJ, Smith DL, Hay SI, Cibulskis RE, Gething PW. The effect of malaria control on \(Plasmodium\, falciparum\) in africa between 2000 and 2015. Nature. 2015;526(7572):207–11.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  5. Davis WA, Clarke PM, Siba PM, Karunajeewa HA, Davy C, Mueller I, Davis TME. Cost-effectiveness of artemisinin combination therapy for uncomplicated malaria in children: data from papua new guinea. Bull World Health Org. 2011;89(3):211–20. https://0-doi-org.brum.beds.ac.uk/10.2471/BLT.10.084103.

    Article  PubMed  PubMed Central  Google Scholar 

  6. Meshnick SR, Taylor TE, Kamchonwongpaisan S. Artemisinin and the antimalarial endoperoxides: from herbal remedy to targeted chemotherapy. Microbiol Rev. 1996;60(2):301–15.

    CAS  PubMed  PubMed Central  Google Scholar 

  7. Shandilya A, Chacko S, Jayaram B, Ghosh I. A plausible mechanism for the antimalarial activity of artemisinin: a computational approach. Sci Rep. 2013;3:2513. https://0-doi-org.brum.beds.ac.uk/10.1038/srep02513.

    Article  PubMed  PubMed Central  Google Scholar 

  8. Sidhu ABS, Verdier-Pinard D, Fidock DA. Chloroquine resistance in \(Plasmodium\, falciparum\) malaria parasites conferred by pfcrt mutations. Science. 2002;298(5591):210–3. https://0-doi-org.brum.beds.ac.uk/10.1126/science.1074045.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  9. Noedl H, Se Y, Schaecher K, Smith BL, Socheat D, Fukuda MM. Evidence of artemisinin-resistant malaria in western cambodia. N Engl J Med. 2008;359(24):2619–20. https://0-doi-org.brum.beds.ac.uk/10.1056/NEJMc0805011.

    Article  CAS  PubMed  Google Scholar 

  10. WHO. Emergency response to artemisinin resistance in the Greater Mekong Subregion. Regional framework for action 2013–2015. Geneva: World Health Organization; 2013. p. 32.

    Google Scholar 

  11. Ataide R, Ashley EA, Powell R, Chan J-A, Malloy MJ, O’Flaherty K, Takashima E, Langer C, Tsuboi T, Dondorp AM, Day NP, Dhorda M, Fairhurst RM, Lim P, Amaratunga C, Pukrittayakamee S, Hien TT, Htut Y, Mayxay M, Faiz MA, Beeson JG, Nosten F, Simpson JA, White NJ, Fowkes FJI. Host immunity to \(Plasmodium\, falciparum\) and the assessment of emerging artemisinin resistance in a multinational cohort. Proc Natl Acad Sci. 2017;114:3515–20. https://0-doi-org.brum.beds.ac.uk/10.1073/pnas.1615875114.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  12. Lubell Y, Dondorp A, Guérin PJ, Drake T, Meek S, Ashley E, Day NP, White NJ, White LJ. Artemisinin resistance—modelling the potential human and economic costs. Malaria J. 2014;13(1):452. https://0-doi-org.brum.beds.ac.uk/10.1186/1475-2875-13-452.

    Article  Google Scholar 

  13. Ashley EA, Dhorda M, Fairhurst RM, Amaratunga C, Lim P, Suon S, Sreng S, Anderson JM, Mao S, Sam B, Sopha C, Chuor CM, Nguon C, Sovannaroth S, Pukrittayakamee S, Jittamala P, Chotivanich K, Chutasmit K, Suchatsoonthorn C, Runcharoen R, Hien TT, Thuy-Nhien NT, Thanh NV, Phu NH, Htut Y, Han K-T, Aye KH, Mokuolu OA, Olaosebikan RR, Folaranmi OO, Mayxay M, Khanthavong M, Hongvanthong B, Newton PN, Onyamboko MA, Fanello CI, Tshefu AK, Mishra N, Valecha N, Phyo AP, Nosten F, Yi P, Tripura R, Borrmann S, Bashraheil M, Peshu J, Faiz MA, Ghose A, Hossain MA, Samad R, Rahman MR, Hasan MM, Islam A, Miotto O, Amato R, MacInnis B, Stalker J, Kwiatkowski DP, Bozdech Z, Jeeyapant A, Cheah PY, Sakulthaew T, Chalk J, Intharabut B, Silamut K, Lee SJ, Vihokhern B, Kunasol C, Imwong M, Tarning J, Taylor WJ, Yeung S, Woodrow CJ, Flegg JA, Das D, Smith J, Venkatesan M, Plowe CV, Stepniewska K, Guerin PJ, Dondorp AM, Day NP, White NJ. Spread of artemisinin resistance in \(Plasmodium\, falciparum\) malaria. N Engl J Med. 2014;371(5):411–23. https://0-doi-org.brum.beds.ac.uk/10.1056/NEJMoa1314981.

    Article  PubMed  PubMed Central  Google Scholar 

  14. WHO/GMP, G.M.P. Insecticide-treated mosquito nets: a WHO position statement. Geneva: World Health Organization; 2011.

  15. Hawley WA, Phillips-Howard PA, ter Kuile FO, Terlouw DJ, Vulule JM, Ombok M, Nahlen BL, Gimnig JE, Kariuki SK, Kolczak MS, Hightower AW. Community-wide effects of permethrin-treated bed nets on child mortality and malaria morbidity in western kenya. Am J Trop Med Hyg. 2003;68(4 Suppl.):121–7.

    PubMed  Google Scholar 

  16. Elissa N, Mouchet J, Riviere F, Meunier JY, Yao K. Resistance of \(Anopheles\, gambiae\) s.s. to pyrethroids in Côte d’Ivoire. Ann Soc Belg Med Trop. 1993;73(4):291–4.

    CAS  PubMed  Google Scholar 

  17. Ranson H, Lissenden N. Insecticide resistance in african anopheles mosquitoes: a worsening situation that needs urgent action to maintain malaria control. Trends Parasitol. 2016;32(3):187–96. https://0-doi-org.brum.beds.ac.uk/10.1016/j.pt.2015.11.010 (Special Issue: Vectors).

    Article  CAS  PubMed  Google Scholar 

  18. Edi CVA, Koudou BG, Jones CM, Weetman D, Ranson H. Multiple-insecticide resistance in \(Anopheles\, gambiae\) mosquitoes, southern Côte d’Ivoire. Emerg Infect Dis. 2012;18(9):1508–11. https://0-doi-org.brum.beds.ac.uk/10.3201/eid1809.120262.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  19. WHO. Pesticide evaluation scheme; 2017. http://www.who.int/whopes/en/. 01 Mar 2017.

  20. Cisse MBM, Keita C, Dicko A, Dengela D, Coleman J, Lucas B, Mihigo J, Sadou A, Belemvire A, George K, Fornadel C, Beach R. Characterizing the insecticide resistance of \(Anopheles\, gambiae\) in Mali. Malaria J. 2015;14(1):327. https://0-doi-org.brum.beds.ac.uk/10.1186/s12936-015-0847-4.

    Article  Google Scholar 

  21. Knipling EF. Possibilities of insect control or eradication through the use of sexually sterile males. J Econ Entomol. 1955;48(4):459. https://0-doi-org.brum.beds.ac.uk/10.1093/jee/48.4.459.

    Article  Google Scholar 

  22. Krafsur E. Sterile insect technique for suppressing and eradicating insect population: 55 years and counting. J Agric Entomol. 1998;15(4):303–17.

    Google Scholar 

  23. Dyck VA, Hendrichs J, Robinson AS. Sterile insect technique: principles and practice in area-wide integrated pest management. Dordrecht: Springer; 2005.

    Book  Google Scholar 

  24. Thistlewood HMA, Judd GJR. Area-wide management of codling moth, Cydia pomonella, in rural and orchard areas. IOBC WPRS Bull. 2003;26(11):103–9.

    Google Scholar 

  25. Thistlewood HMA, Judd GJR, Clodius MEO. Sustainable management and monitoring of codling moth, \(Cydia\, pomonella\) in an area-wide program employing sit. In: IAEA, editor. Improvement of codling moth SIT to facilitate expansion of field application. Working Paper series. Vienna: Atomic Energy Agency; 2004. p. 79– 87.

  26. Robinson AS, Hendrichs J. Prospects for the future development and application of the sterile insect technique. In: Dyck VA, Hendrichs J, Robinson AS, editors. Sterile insect technique: principles and practice in area-wide integrated pest management. Berlin: Springer; 2005. p. 727–60.

    Chapter  Google Scholar 

  27. Esteva L, Yang HM. Mathematical model to assess the control of \(Aedes\, aegypti\) mosquitoes by the sterile insect technique. Math Biosci. 2005;198(2):132–47.

    Article  PubMed  Google Scholar 

  28. White SM, Rohani P, Sait SM. Modelling pulsed releases for sterile insect techniques: fitness costs of sterile and transgenic males and the effects on mosquito dynamics. J Appl Ecol. 2010;47(6):1329–39.

    Article  Google Scholar 

  29. Stone CM. Transient population dynamics of mosquitoes during sterile male releases: modelling mating behaviour and perturbations of life history parameters. PLoS ONE. 2013;8(9):1–13. https://0-doi-org.brum.beds.ac.uk/10.1371/journal.pone.0076228.

    Article  Google Scholar 

  30. Seirin Lee S, Baker RE, Gaffney EA, White SM. Modelling \(Aedes\, aegypti\) mosquito control via transgenic and sterile insect techniques: endemics and emerging outbreaks. J Theor Biol. 2013;331:78–90. https://0-doi-org.brum.beds.ac.uk/10.1016/j.jtbi.2013.04.014.

    Article  CAS  PubMed  Google Scholar 

  31. Muriu SM, Coulson T, Mbogo CM, Godfray HCJ. Larval density dependence in \(Anopheles\, gambiae\) s.s., the major african vector of malaria. J Animal Ecol. 2013;82(1):166–74. https://0-doi-org.brum.beds.ac.uk/10.1111/1365-2656.12002.

    Article  Google Scholar 

  32. Phuc HK, Andreasen MH, Burton RS, Vass C, Epton MJ, Pape G, Fu G, Condon KC, Scaife S, Donnelly CA, Coleman PG, White-Cooper H, Alphey L. Late-acting dominant lethal genetic systems and mosquito control. BMC Biol. 2007;5(1):11. https://0-doi-org.brum.beds.ac.uk/10.1186/1741-7007-5-11.

    Article  PubMed  PubMed Central  Google Scholar 

  33. Yakob L, Alphey L, Bonsall MB. \(Aedes\, aegypti\) control: the concomitant role of competition, space and transgenic technologies. J Appl Ecol. 2008;45(4):1258–65. https://0-doi-org.brum.beds.ac.uk/10.1111/j.1365-2664.2008.01498.x.

    Article  Google Scholar 

  34. Alphey N, Bonsall MB. Interplay of population genetics and dynamics in the genetic control of mosquitoes. J R Soc Interface. 2014;11(93):20131071. https://0-doi-org.brum.beds.ac.uk/10.1098/rsif.2013.1071.

    Article  PubMed  PubMed Central  Google Scholar 

  35. Thomas DD, Donnelly CA, Wood RJ, Alphey LS. Insect population control using a dominant, repressible, lethal genetic system. Science. 2000;287(5462):2474–6.

    Article  CAS  PubMed  Google Scholar 

  36. Alphey L, Nimmo D, O’Connell S, Alphey N. Insect population suppression using engineered insects. In: Aksoy S, editor. Transgenesis and the management of vector-borne disease. New York: Springer; 2008. p. 93–103.

    Chapter  Google Scholar 

  37. Atkinson MP, Su Z, Alphey N, Alphey LS, Coleman PG, Wein LM. Analyzing the control of mosquito-borne diseases by a dominant lethal genetic system. Proc Natl Acad Sci USA. 2007;104(22):9540–5. https://0-doi-org.brum.beds.ac.uk/10.1073/pnas.0610685104.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  38. Sinclair D, Zani B, Donegan S, Olliaro P, Garner P. Artemisinin-based combination therapy for treating uncomplicated malaria. Cochrane Database Syst Rev. 2009;. Art. No. CD007483. https://0-doi-org.brum.beds.ac.uk/10.1002/14651858.CD007483.pub2.

  39. Lengeler C. Insecticide-treated bednets and curtains for preventing malaria. Cochrane Database Syst Rev. 2004; Art. No. CD000363. https://0-doi-org.brum.beds.ac.uk/10.1002/14651858.CD000363.pub2.

  40. Pluess B, Tanser FC, Lengeler C, Sharp BL. Indoor residual spraying for preventing malaria. Cochrane Database Syst Rev. 2010;. Art. No.: CD006657. https://0-doi-org.brum.beds.ac.uk/10.1002/14651858.CD006657.pub2.

  41. Goodman CA, Coleman PG, Mills AJ. Cost-effectiveness of malaria control in sub-saharan africa. Lancet. 1999;354(9176):378–85. https://0-doi-org.brum.beds.ac.uk/10.1016/S0140-6736(99)02141-8.

    Article  CAS  PubMed  Google Scholar 

  42. Mumford JD. Application of benefit/cost analysis to insect pest control using the sterile insect technique. In: Dyck VA, Hendrichs J, Robinson AS, editors. Sterile insect technique: principles and practice in area-wide integrated pest management. Berlin: Springer; 2005. p. 481–98.

    Chapter  Google Scholar 

  43. Alphey N, Alphey L, Bonsall MB. A model framework to estimate impact and cost of genetics-based sterile insect methods for dengue vector control. PLoS ONE. 2011;6(10):1–12. https://0-doi-org.brum.beds.ac.uk/10.1371/journal.pone.0025384.

    Article  Google Scholar 

  44. Gurney WSC, Blythe SP, Nisbet RM. Nicholson’s blowflies revisited. Nature. 1980;287:17–21.

    Article  Google Scholar 

  45. Gurney WSC, Nisbet RM, Lawton JH. The systematic formulation of tractable single-species population models incorporating age structure. J Animal Ecol. 1983;52(2):479–95.

    Article  Google Scholar 

  46. Ross R. The prevention of malaria. London: J. Murray; 1910.

    Google Scholar 

  47. Ross R. An application of the theory of probabilities to the study of a priori pathometry. Part i. Proc R Soc London A Math Phys Eng Sci. 1916;92(638):204–30. https://0-doi-org.brum.beds.ac.uk/10.1098/rspa.1916.0007.

    Article  Google Scholar 

  48. MacDonald G. The analysis of equilibrium in malaria. Trop Dis Bull. 1952;49(9):813–29.

    CAS  PubMed  Google Scholar 

  49. MacDonald G. The Epidemiology and control of malaria. Oxford: Oxford Medical Publications. Oxford University Press; 1957.

    Google Scholar 

  50. Aron JL, May RM. The population dynamics of malaria. In: Anderson RM, editor. Population dynamics of infectious diseases: theory and application. Population and community biology. London: Chapman and Hall; 1982. p. 139–79.

    Chapter  Google Scholar 

  51. WHO. Global health observatory data repository. Geneva: World Health Organization; 2018. http://apps.who.int/gho/data/. 01 Mar 2017.

  52. Kilama M, Smith DL, Hutchinson R, Kigozi R, Yeka A, Lavoy G, Kamya MR, Staedke SG, Donnelly MJ, Drakeley C, Greenhouse B, Dorsey G, Lindsay SW. Estimating the annual entomological inoculation rate for \(Plasmodium\, falciparum\) transmitted by \(Anopheles\, gambiae\) s.l. using three sampling methods in three sites in uganda. Malaria J. 2014;13(1):111. https://0-doi-org.brum.beds.ac.uk/10.1186/1475-2875-13-111.

    Article  Google Scholar 

  53. Takken W, Klowden MJ, Chambers GM. Effect of body size on host seeking and blood meal utilization in \(Anopheles\, gambiae\) sensu stricto (diptera: Culicidae): the disadvantage of being small. J Med Entomol. 1998;35(5):639–45.

    Article  CAS  PubMed  Google Scholar 

  54. Bellows TS. The descriptive properties of some models for density dependence. J Animal Ecol. 1981;50(1):139–56.

    Article  Google Scholar 

  55. Maynard Smith J, Slatkin M. The stability of predator-prey systems. Ecology. 1973;54(2):384–91. https://0-doi-org.brum.beds.ac.uk/10.2307/1934346.

    Article  Google Scholar 

  56. Torgerson DJ, Raftery J. Discounting. BMJ. 1999;319(7214):914–5. https://0-doi-org.brum.beds.ac.uk/10.1136/bmj.319.7214.914.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  57. Clark CW. Mathematical bioeconomics: the mathematics of conservation. Pure and applied mathematics: a Wiley series of texts, monographs and tracts. Hoboken: John Wiley & Sons; 2010.

    Google Scholar 

  58. Hurford A, Cownden D, Day T. Next-generation tools for evolutionary invasion analyses. J R Soc Interface. 2010;7(45):561–71. https://0-doi-org.brum.beds.ac.uk/10.1098/rsif.2009.0448.

    Article  PubMed  Google Scholar 

  59. May RM, Anderson RM. Endemic infections in growing populations. Math Biosci. 1985;77:141–56.

    Article  Google Scholar 

  60. Diekmann O, Heesterbeek JAP, Metz JAJ. On the definition and the computation of the basic reproduction ratio \(r_0\) in models for infectious diseases in heterogeneous populations. J Math Biol. 1990;28(4):365–82.

    Article  CAS  PubMed  Google Scholar 

  61. White MT, Conteh L, Cibulskis R, Ghani AC. Costs and cost-effectiveness of malaria control interventions—a systematic review. Malaria J. 2011;10(1):337. https://0-doi-org.brum.beds.ac.uk/10.1186/1475-2875-10-337.

    Article  Google Scholar 

  62. Khamis SH. A new system of index numbers for national and international purposes. J R Stat Soc Ser A Gen. 1972;135(1):96–121.

    Article  Google Scholar 

  63. Taylor C, Touré YT, Carnahan J, Norris DE, Dolo G, Traoré SF, Edillo FE, Lanzaro GC. Gene flow among populations of the malaria vector, \(Anopheles\, gambiae\), in Mali, West Africa. Genetics. 2001;157(2):743–50.

    CAS  PubMed  PubMed Central  Google Scholar 

  64. Hendron R-WS, Bonsall MB. The interplay of vaccination and vector control on small dengue networks. J Theor Biol. 2016;407:349–61. https://0-doi-org.brum.beds.ac.uk/10.1016/j.jtbi.2016.07.034.

    Article  PubMed  PubMed Central  Google Scholar 

  65. Okell LC, Drakeley CJ, Bousema T, Whitty CJM, Ghani AC. Modelling the impact of artemisinin combination therapy and long-acting treatments on malaria transmission intensity. PLoS Med. 2008;5(11):1–12. https://0-doi-org.brum.beds.ac.uk/10.1371/journal.pmed.0050226.

    Article  Google Scholar 

  66. Hastings IM, Watkins WM, White NJ. The evolution of drug-resistant malaria: the role of drug elimination half-life. Philos Trans R Soc B Biol Sci. 2002;357(1420):505–19. https://0-doi-org.brum.beds.ac.uk/10.1098/rstb.2001.1036.

    Article  CAS  Google Scholar 

  67. Hastings IM, Hodel EM. Pharmacological considerations in the design of anti-malarial drug combination therapies—is matching half-lives enough? Malaria J. 2014;13(1):62. https://0-doi-org.brum.beds.ac.uk/10.1186/1475-2875-13-62.

    Article  Google Scholar 

  68. Tyson R, Newton KD, Thistlewood HMA, Judd GJR. Mating rates between sterile and wild codling moths (\(Cydia\, pomonella\)) in springtime: a simulation study. J Theor Biol. 2008;254:319–30.

    Article  CAS  PubMed  Google Scholar 

  69. Alegana VA, Kigozi SP, Nankabirwa J, Arinaitwe E, Kigozi R, Mawejje H, Kilama M, Ruktanonchai NW, Ruktanonchai CW, Drakeley C, Lindsay SW, Greenhouse B, Kamya MR, Smith DL, Atkinson PM, Dorsey G, Tatem AJ. Spatio-temporal analysis of malaria vector density from baseline through intervention in a high transmission setting. Parasites Vectors. 2016;9(1):637. https://0-doi-org.brum.beds.ac.uk/10.1186/s13071-016-1917-3.

    Article  PubMed  PubMed Central  Google Scholar 

  70. Bonsall MB, Yakob L, Alphey N, Alphey L. Transgenic control of vectors: the effects of interspecific interactions. Isr J Ecol Evol. 2010;56:353–70.

    Article  Google Scholar 

  71. Trape J-F, Lefebvre-Zante E, Legros F, Ndiaye G, Bouganali H, Druilhe P, Salem G. Vector density gradients and the epidemiology of urban malaria in Dakar, Senegal. Am J Trop Med Hyg. 1992;47(2):181–9. https://0-doi-org.brum.beds.ac.uk/10.4269/ajtmh.1992.47.181.

    Article  CAS  PubMed  Google Scholar 

  72. Ototo EN, Mbugi JP, Wanjala CL, Zhou G, Githeko AK, Yan G. Surveillance of malaria vector population density and biting behaviour in Western Kenya. Malaria J. 2015;14(1):244. https://0-doi-org.brum.beds.ac.uk/10.1186/s12936-015-0763-7.

    Article  Google Scholar 

  73. Molineaux L, Gramiccia G. The garki project: research on the epidemiology and control of malaria in the Sudan Savanna of West Africa. Geneva: World Health Organization; 1980.

    Google Scholar 

  74. Newton EA, Reiter P. A model of the transmission of dengue fever with an evaluation of the impact of ultra-low volume (ULV) insecticide applications on dengue epidemics. Am J Trop Med Hyg. 1992;47(6):709–20.

    Article  CAS  PubMed  Google Scholar 

  75. Churcher TS, Sinden RE, Edwards NJ, Poulton ID, Rampling TW, Brock PM, Griffin JT, Upton LM, Zakutansky SE, Sala KA, Angrisano F, Hill AVS, Blagborough AM. Probability of transmission of malaria from mosquito to human is regulated by mosquito parasite density in naïve and vaccinated hosts. PLoS Pathog. 2017;13(1):1–18. https://0-doi-org.brum.beds.ac.uk/10.1371/journal.ppat.1006108.

    Article  Google Scholar 

  76. Ermert V, Fink AH, Jones AE, Morse AP. Development of a new version of the liverpool malaria model. i. refining the parameter settings and mathematical formulation of basic processes based on a literature review. Malaria J. 2011;10(1):35. https://0-doi-org.brum.beds.ac.uk/10.1186/1475-2875-10-35.

    Article  Google Scholar 

  77. Ermert V, Fink AH, Jones AE, Morse AP. Development of a new version of the liverpool malaria model. II. Calibration and validation for West Africa. Malaria J. 2011;10(1):62. https://0-doi-org.brum.beds.ac.uk/10.1186/1475-2875-10-62.

    Article  Google Scholar 

  78. Ajayi IO, Browne EN, Bateganya F, Yar D, Happi C, Falade CO, Gbotosho GO, Yusuf B, Boateng S, Mugittu K, Cousens S, Nanyunja M. Effectiveness of artemisinin-based combination therapy used in the context of home management of malaria: a report from three study sites in sub-saharan africa. Malaria J. 2008;7:190.

    Article  Google Scholar 

  79. Dye C. Models for the population dynamics of the yellow fever mosquito, \(Aedes\, aegypti\). J Animal Ecol. 1984;53(1):247–68.

    Article  Google Scholar 

  80. Southwood TRE, Murdie G, Yasuno M, Tonn RJ, Reader PM. Studies on the life budget of \(Aedes\, aegypti\) in Wat Samphaya, Bangkok, Thailand. Bull World Health Org. 1972;46(2):211–26.

    CAS  PubMed  PubMed Central  Google Scholar 

  81. Sheppard PM, Macdonald WW, Tonn RJ, Grab B. The dynamics of an adult population of \(Aedes\, aegypti\) in relation to dengue haemorrhagic fever in bangkok. J Animal Ecol. 1969;38(3):661–702.

    Article  Google Scholar 

  82. Hoshen MB, Morse AP. A weather-driven model of malaria transmission. Malaria J. 2004;3(1):32. https://0-doi-org.brum.beds.ac.uk/10.1186/1475-2875-3-32.

    Article  Google Scholar 

  83. Deredec A, Godfray HCJ, Burt A. Requirements for effective malaria control with homing endonuclease genes. Proc Natl Acad Sci. 2011;108(43):874–80. https://0-doi-org.brum.beds.ac.uk/10.1073/pnas.1110717108.

    Article  Google Scholar 

  84. Clements AN, Paterson GD. The analysis of mortality and survival rates in wild populations of mosquitoes. J Appl Ecol. 1981;18(2):373–99.

    Article  Google Scholar 

  85. Jannat KN-E, Roitberg BD. Effects of larval density and feeding rates on larval life history traits in \(Anopheles\, gambiae\) s.s. (diptera: Culicidae). J Vector Ecol. 2013;38(1):120–6.

    Article  PubMed  Google Scholar 

  86. Parker A. Mass-rearing for sterile insect release. In: Dyck VA, Hendrichs J, Robinson AS, editors. Sterile insect technique: principles and practice in area-wide integrated pest management. Dordrecht: Springer; 2005. p. 209–32.

    Chapter  Google Scholar 

  87. Bailey DL, Lowe RE, Dame DA, Seawright JA. Mass rearing the genetically altered macho strain of \(Anopheles\, albimanus\) Wiedemann. Am J Trop Med Hyg. 1980;29(1):141–9.

    Article  CAS  PubMed  Google Scholar 

  88. Dame DA. Genetic control by sterilized mosquitoes. In: Chapman HCEA, editor. Biological control of mosquitoes. Township: American Mosquito Control Association; 1985. p. 159–72.

    Google Scholar 

  89. Singh KRP, Razdan RK. Mass rearing of Culex pipiens fatigans Wied. Under ambient conditions. WHO VBC. 1975;75:537.

    Google Scholar 

  90. International Atomic Energy Agency. Model Business Plan for a Sterile Insect Production Facility. Vienna: IAEA; 2008.

    Google Scholar 

  91. Lenhart S, Workman JT. Optimal control applied to biological models. Chapman & Hall/CRC mathematical and computational biology. Abingdon: Taylor & Francis; 2007.

    Google Scholar 

Download references

Authors’ contributions

DK, CEM, KK and MBB conceived the study. DK performed the mathematical analysis and drafted the manuscript. CEM, KK and MBB helped draft the manuscript. All authors gave final approval for publication. All authors read and approved the final manuscript.

Acknowledgements

We thank an anonymous reviewer, Luke Alphey, Jack Newman and the Mathematical Ecology Research Group for their helpful comments and advice on this work.

Competing interests

The authors declare that they have no competing interests.

Availability of data and materials

Data and code are available through the Open Science Framework at http://osf.io/8kf6x.

Consent for publication

Not applicable.

Ethics approval and consent to participate

Not applicable.

Funding

The authors are grateful to DARPA for financial support (Contract: HR00111-16-2-0005).

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Author information

Authors and Affiliations

Authors

Corresponding author

Correspondence to Michael B. Bonsall.

Appendices

Appendix A: Optimal control formulation

Optimal control strategies for disease control rely on defining objective functions from which the cost-effectiveness of interventions can be evaluated. Here the costs of treating infected individuals, the costs of drug based therapies and the costs of SIT mosquito releases are explored, and an objective functional is defined to explore different cost structures (costs are linear, costs are convex, costs are concave).

Two control parameters are present in the model (1): the release ratio u and the proportion medicated w as optimal control strategy to manage the spread of malaria is required. This strategy would minimize the socioeconomic burden of malaria, the cost of artemisinin medication and the cost of vector control via sterile insect release. If the cost of direct and indirect health care and lost productivity due to a single case of malaria is \(\theta _{1}\) per day, then the daily cost of having a proportion h of a population of size N infected is \(\theta _{1} N h,\) which has units of dollars per day. To capture the potential nonlinearity of cost, the cost function is defined as \(C_{h} = \theta _{1} N h^{m_{1}},\) which also has units of dollars per day. The exponent \(0<m_{1}\le 2\) affects only the dimensionless quantity h, and varying \(m_1\) over its range allows us to investigate cost structures that are linear, convex or concave. Similarly, if the dollar cost of medicating a single patient with artemisinin-based therapy is \(\theta _{2}\) per day, then the daily cost of medicating a proportion w of an infected population of size Nh is \(\theta _{2} N h w.\) As before, a cost function is defined, \(C_{w} = \theta _{2} N h w^{m_{2}},\) which has units of dollars per day and where the exponent \(0<m_{2}\le 2\) affects only the dimensionless quantity w. Finally, if the cost of rearing and releasing one (or 100, or 1000, depending on population units) sterile or GM insect per day is \(\theta _{3},\) then releasing a proportion u of the equilibrium vector population \(A^{*}\) daily will cost \(\theta _{3} A^{*} u\) per day. Again, to allow for possible nonlinearity in costing, the cost function is defined as \(C_{u} = \theta _{3} A^{*} u^{m_{3}},\) which has units of dollars per day independent of the value of the exponent \(0<m_{3}\le 2,\) which affects only the dimensionless quantity u.

Using the cost functions defined above, a cost functional is defined

$$\begin{aligned} J[\varvec{x}(u,w)] = \int _0^T \mathrm {e}^{-\psi t} \left[ \theta _{1} N h^{m_{1}} + \theta _{2} N h w^{m_{2}} + \theta _{3} A^{*} u^{m_{3}}\right] \mathrm {d}t, \end{aligned}$$
(8)

where future costs are discounted at a constant rate \(\psi\) to deal with opportunity costs in investing in public health control for immediate benefits [56], and \(\varvec{x} = (L, A, h, v)^{\top }\) is the vector of state variables. The cost function J measures the total cost of disease, artemisinin medication and vector control over a fixed time period T, and has units of dollars. From (8) a new functional is defined \(J_{pc} = J/N\), which measures the total cost per capita:

$$\begin{aligned} J_{pc}[\varvec{x}(u,w)] = \int _0^T \mathrm {e}^{-\psi t} \left[ \theta _{1} h^{m_{1}} + \theta _{2} h w^{m_{2}} + \theta _{3} k^{*} u^{m_{3}}\right] \mathrm {d}t, \end{aligned}$$
(9)

where \(k^{*} = A^{*}/N\) is the number of vectors per host at control-free equilibrium. The goal, then, is to minimize the total cost of the malaria management strategy:

$$\begin{aligned} \mathop {\mathop {\min }\limits _{u \ge 0} }\limits _{0 \le w \le 0.6} \{ {J_{pc}}\} . \end{aligned}$$
(10)

To solve the optimization problem (1) (9) and (10), the Hamiltonian and adjoint system (see below) are formed and analytical expressions for the characterizations of the optimal controls \(u^*(t)\) and \(w^*(t)\) are derived that minimize the Hamiltonian (see Appendix A.1: “Optimality conditions”). To solve the system, the forward–backward sweep method is used with Runge–Kutta integration schemes [91].

For early-acting SIT, the Hamiltonian is

$$\begin{aligned} H =&\,\, \mathrm {e}^{-\psi t} \left[ \theta _{1} h^{m_{1}} + \theta _{2} h w^{m_{2}} + \theta _{3} k^{*} u^{m_{3}}\right] + \lambda _{1}\left[ \rho A C(A,u)\nonumber \right. \\&\left. - f(L) L - (m + \mu _L) L\right] + \lambda _{2}\left[ \frac{m}{2} L - \mu _A A\right] \nonumber \\&\,+ \lambda _{3}\left[ \frac{A}{N} b a (1 - h) v - (\gamma + w s) h\right] + \lambda _{4}\left[ b c (1 - v) h - \frac{m L}{2 A} v\right] , \end{aligned}$$
(11)

where \(\lambda _{i}(t)\) are the adjoint variables. Each adjoint variable \(\lambda _{i}\) satisfies an equation found by differentiating the Hamiltonian with respect to its corresponding state variable \(x_{i}\) and negating:

$$\begin{aligned} \frac{\mathrm {d}\lambda _1}{\mathrm {d}t} =&\, -\frac{\partial H}{\partial L} = \lambda _{1}(f'(L)L + f(L) + m + \mu _{L}) - \frac{1}{2}\lambda _{2}m + \lambda _{4}\frac{m v}{2 A}, \end{aligned}$$
(12a)
$$\begin{aligned} \frac{\mathrm {d}\lambda _{2}}{\mathrm {d}t} =&\, -\frac{\partial H}{\partial A} = -\lambda _{1}\rho C(A,u) \left( 1 + A \frac{\partial }{\partial A}\log {C(A,u)}\right) \nonumber \\&+ \lambda _{2}\mu _{A} - \lambda _{3}\frac{b a}{N}(1 - h) v - \lambda _{4}\frac{m L v}{2 A^{2}}, \end{aligned}$$
(12b)
$$\begin{aligned} \frac{\mathrm {d}\lambda _{3}}{\mathrm {d}t} =&\, -\frac{\partial H}{\partial h} = -\mathrm {e}^{-\psi t}(m_{1} \theta _{1} h^{m_{1}-1} + \theta _{2} w^{m_{2}}) + \lambda _{3}\left( \frac{A}{N} b a v + \gamma + w s\right) - \lambda _{4}b c (1 - v), \end{aligned}$$
(12c)
$$\begin{aligned} \frac{\mathrm {d}\lambda _{4}}{\mathrm {d}t} =&\, -\frac{\partial H}{\partial v} = -\lambda _{3}\frac{A}{N} b a (1 - h) + \lambda _{4}\left( b c h + \frac{m L}{2 A}\right) . \end{aligned}$$
(12d)

The system (12) is closed by the transversality condition

$$\begin{aligned} \lambda _{i}(T) = 0. \end{aligned}$$
(13)

For the late-acting SIT, which prevents larvae from reaching the adult stage [the dynamics of which are captured in (5)], the Hamiltonian is

$$\begin{aligned} H =&\,\, \mathrm {e}^{-\psi t} \left[ \theta _{1} h^{m_{1}} + \theta _{2} h w^{m_{2}} + \theta _{3} k^{*} u^{m_{3}}\right] + \lambda _{1}\left[ \rho A - f(L) L - (m + \mu _L) L\right] \nonumber \\&+\lambda _{2}\left[ \frac{m}{2} L C(A,u) - \mu _A A\right] \nonumber \\&\,+ \lambda _{3}\left[ \frac{A}{N} b a (1 - h) v - (\gamma + w s) h\right] + \lambda _{4}\left[ b c (1 - v) h - \frac{m L}{2 A} C(A,u) v\right] , \end{aligned}$$
(14)

and the adjoint system is

$$\begin{aligned} \frac{\mathrm {d}\lambda _1}{\mathrm {d}t} =&\, -\frac{\partial H}{\partial L} = \lambda _{1}(f'(L)L + f(L) + m + \mu _{L}) - \frac{1}{2} \lambda _{2} m C(A,u) + \lambda _{4}\frac{m}{2 A} C(A,u) v, \end{aligned}$$
(15a)
$$\begin{aligned} \frac{\mathrm {d}\lambda _{2}}{\mathrm {d}t} &=-\frac{\partial H}{\partial A} = -\lambda _{1}\rho - \lambda _{2}\left( \frac{m}{2} L \frac{\partial }{\partial A}C(A,u) - \mu _{A}\right) - \lambda _{3}\frac{b a}{N}(1 - h) v \nonumber \\& \quad - \lambda _{4}\frac{m L v}{2 A^{2}}C(A,u)\left( 1 - A \frac{\partial }{\partial A}\log {C(A,u)}\right) , \end{aligned}$$
(15b)
$$\begin{aligned} \frac{\mathrm {d}\lambda _{3}}{\mathrm {d}t} =&\, -\frac{\partial H}{\partial h} = -\mathrm {e}^{-\psi t}(m_{1} \theta _{1} h^{m_{1}-1} + \theta _{2} w^{m_{2}}) + \lambda _{3}\left( \frac{A}{N} b a v + \gamma + w s\right) - \lambda _{4}b c (1 - v), \end{aligned}$$
(15c)
$$\begin{aligned} \frac{\mathrm {d}\lambda _{4}}{\mathrm {d}t} =&\, -\frac{\partial H}{\partial v} = -\lambda _{3}\frac{A}{N} b a (1 - h) + \lambda _{4}\left( b c h + \frac{m L}{2 A}C(A,u)\right) . \end{aligned}$$
(15d)

with boundary conditions as defined in (13).

A.1 Optimality conditions

The optimal controls \((u^{*}, w^{*})\) minimize the Hamiltonian H. Thus,

$$\begin{aligned} \frac{\partial H}{\partial u} = 0 \quad \mathrm {at} \,\, u = u^{*},\quad\frac{\partial H}{\partial w} = 0 \quad \mathrm {at} \,\, w = w^{*}. \end{aligned}$$
(16)

For the early-acting SIT, artemisinin control is investigated. When \(m_{2}\ne 1,\)

$$\begin{aligned} \frac{\partial H}{\partial w} = m_{2} \theta _{2} h w^{m_{2}-1}\mathrm {e}^{-\psi t} - \lambda _{3} s h = 0 \qquad \mathrm {at} \,\, w = w^{*}. \end{aligned}$$
(17)

The optimal control \(w^{*}\) may then be found:

$$\begin{aligned} w^{*} = \left( \frac{\lambda _{3} s}{m_{2} \theta _{2}}\mathrm {e}^{\psi t}\right) ^{\frac{1}{m_{2}-1}}. \end{aligned}$$
(18)

The second derivative of H indicates whether the solution is a maximum or a minimum:

$$\begin{aligned} \frac{\partial ^2 H}{\partial w ^2} = (m_2 - 1) m_2 \theta _2 h w^{m_2 - 2} \mathrm {e}^{\psi t}; \end{aligned}$$
(19)

this is only positive if \(m_2>1,\) so the solution (18) is restricted to cost functions with \(m_2>1.\)

For the vector control, consider the case \(m_{3}=1.\) Differentiating the Hamiltonian for early-acting SIT produces

$$\begin{aligned} \frac{\partial H}{\partial u} = m_{3} \theta _{3}\frac{A^{*}}{N}\mathrm {e}^{-\psi t} - \frac{\lambda _{1} \rho A^{2} A^{*}}{(A + A^{*}u)^{2}} = 0 \qquad \mathrm {at} \,\, u = u^{*}, \end{aligned}$$
(20)

which may be rearranged to form a quadratic in u and subsequently solved to give

$$\begin{aligned} u^{*} = - A_{r} \pm \sqrt{\chi }, \end{aligned}$$
(21)

where \(A_{r} = A/A^{*}\) and

$$\begin{aligned} \chi = A_{r}^{2} \frac{\lambda _{1} \rho N}{m_{3} \theta _{3}}\mathrm {e}^{\psi t}. \end{aligned}$$
(22)

The concavity condition holds for (21), as \(\frac{\partial ^2 H}{\partial u ^2}>0.\) For late-acting SIT, the form (21) still holds but with

$$\begin{aligned} \chi = \frac{m L A_{r}^{2}A^{*}}{2 m_{3} \theta _{3} k^{*}}\mathrm {e}^{\psi t} ( \lambda _2 A - \lambda _4 v ). \end{aligned}$$
(23)

For the case \(m_{3}\ne 1,\) differentiating the Hamiltonian leads to a rather more complicated equation for the characterization of \(u^{*.}\) For the quadratic cost function \(m_3 = 2,\) the derivative of the Hamiltonian for early-acting SIT is

$$\begin{aligned} \frac{\partial H}{\partial u} = 2\theta _{3} u \frac{A^{*}}{N}\mathrm {e}^{-\psi t} - \frac{\lambda _{1} \rho A^{2} A^{*}}{(A + A^{*}u)^{2}} = 0 \qquad \mathrm {at} \,\, u = u^{*}. \end{aligned}$$
(24)

Equation (24) may be rearranged to form a cubic in u:

$$\begin{aligned} u^{3} + 2 A_{r} u^{2} + A_{r}^{2} u - \chi = 0 \qquad \mathrm {at} \,\, u = u^{*}. \end{aligned}$$
(25)

To aid analysis, (25) is rewritten as a depressed cubic by making the change of variables \(\hat{u} = u + 2A_{r}/3;\) then,

$$\begin{aligned} \hat{u}^{3} + p \hat{u} + q = 0 \qquad \mathrm {at} \;\; u = u^{*}, \end{aligned}$$
(26)

where

$$\begin{aligned} p = -\frac{A_{r}^{2}}{3},\quad q = -\frac{2}{27} A_{r}^{3} - \chi , \end{aligned}$$
(27)

where \(\chi\) is as defined in (22). The three roots of the equation (26) may be succinctly expressed in the form

$$\begin{aligned} \hat{u}^{*}_{j} = 2\sqrt{\frac{-p}{3}}\cos {\left( \frac{1}{3}\arccos {\left( \frac{3 q}{2 p}\sqrt{\frac{3}{-p}}\right) } - \frac{2 \pi j}{3}\right) }, \quad j = 0,1,2, \end{aligned}$$
(28)

from which the optimal control \(u^{*}\) may be constructed by inverting the change of variables,

$$\begin{aligned} u^{*} = \hat{u}^{*} - \frac{2}{3}A_{r}. \end{aligned}$$
(29)

For late-acting SIT the characterization is also defined by (29) with (27) and (28), but with \(\chi\) now defined by (23).

The solutions (29) and (28) will either be all real, or be made up of one real root and two complex conjugate roots; on top of this, any real root may be negative. The character of each root—each value of j in (28)—is not fixed in time, and each may shift from real to complex as the parameters p, q evolve. When choosing the optimal solution, at each point in time, the non-zero, non-negative real root is selected. If (all) the real root(s) is (are) negative at a given time, the control is set to zero at that time. If all three roots are real and positive, the minimum value is taken. For \(m_3 = 2,\) \(\frac{\partial ^2 H}{\partial u ^2}>0\) is found so the solution is a minimum.

An analytical solution may also be found for a square root cost function, \(m_3=\frac{1}{2}.\) In this case, for early-acting SIT, the solution is

$$\begin{aligned} \frac{\partial H}{\partial u} = \frac{1}{2}\theta _{3} u^{-\frac{1}{2}} k^{*}\mathrm {e}^{-\psi t} - \frac{\lambda _{1} \rho A^{2}}{A^{*}(A_r + u)^{2}} = 0 \qquad \mathrm {at} \,\, u = u^{*}. \end{aligned}$$
(30)

Making the change of variables \(\mu = u^{\frac{1}{2}}\) (30) can be written as a depressed quartic equation for \(\mu,\)

$$\begin{aligned} \mu ^{4} + p \mu ^{2} + q \mu + r = 0, \end{aligned}$$
(31)

where

$$\begin{aligned} p = 2 A_r, \quad q = -2 A_{r}^{2}\chi , \quad r = A_{r}^{2}, \end{aligned}$$
(32)

where \(\chi\) is as defined in (22). The solutions may then be written down using the canonical theory,

$$\begin{aligned} \mu = \pm _{1} \frac{1}{2}\sqrt{2 \bar{m}} \pm _{2} \sqrt{-\left( \frac{p + \bar{m}}{2} \pm _{1} \frac{\sqrt{2} q}{4 \sqrt{\bar{m}}}\right) }, \end{aligned}$$
(33)

where

$$\begin{aligned} \bar{m} = \frac{1}{2} \left( -\frac{2}{3} p + \frac{1}{3}\left( Q + \frac{\Delta _0}{Q} \right) \right) , \end{aligned}$$
(34)

with

$$\begin{aligned} Q = \left[ \frac{1}{2}\left( \Delta _1 + \sqrt{\Delta _1^{2} - 4 \Delta _0^{3}}\right) \right] ^{\frac{1}{3}}, \end{aligned}$$
(35)

and

$$\begin{aligned} \Delta _0 = p^2 + 12 r,\quad\Delta _1 = 2 p^3 + 27 q^2 - 72 p r. \end{aligned}$$
(36)

The characterization of the control may then be found by reversing the change of variables, \(u^{*} = \mu ^{2}\). For late-acting SIT, the solution may also be expressed in the form (33), but with \(\chi\) as defined in (23).

In the \(m_{3} = \frac{1}{2}\) case, the concavity condition on \(\partial ^2 H / \partial u^2\) holds only partially throughout the domain. It is found that

$$\begin{aligned} \frac{\partial ^2 H}{\partial u ^2} = \frac{2 \lambda _1 \rho A^2 A^{*^{2}}}{(A + u A^*)^3} - \frac{1}{4}\frac{\theta _3 k^*}{u^{\frac{3}{2}}}\mathrm {e}^{-\psi t}, \end{aligned}$$
(37)

which is positive for large enough u, but goes negative when

$$\begin{aligned} u^{\frac{3}{2}} \lesssim \frac{1}{8}\frac{\theta _3 A}{\lambda _1 \rho A^* N}\mathrm {e}^{-\psi t}. \end{aligned}$$
(38)

This relation generally corresponds to a very low value of control. Since small errors in a near-zero release ratio do not alter the overall cost or efficacy of the control strategy to an appreciable degree, this problem may be ignored. (A similar result is found in the case of late-acting SIT).

Appendix B: \(R_0\) and equilibria

B.1 Basic reproductive ratio

For a constant mosquito population, the next generation matrix is given (using the \(A = F - V,\) \(G=FV^{-1}\) decomposition method of [58]) by

$$\begin{aligned} G = \left( \begin{matrix}0 &{} \frac{k^* b a}{\mu _A} \\ \frac{b c}{\gamma + w s} &{} 0 \end{matrix}\right) , \end{aligned}$$
(39)

for the rare initial infection scenario \((1-h)\sim 1,\) \((1-v)\sim 1.\) The basic reproductive ratio of disease (host to host) is given by the (square of the) dominant eigenvalue of G;

$$\begin{aligned} \begin{vmatrix} -\lambda&\frac{k^* b a}{\mu _A} \\ \frac{b c}{\gamma + w s}&-\lambda \end{vmatrix} = \lambda ^{2} - \frac{k^* b^2 a c}{\mu _A (\gamma + w s)} = 0, \end{aligned}$$
(40)

hence

$$\begin{aligned} R_0 = \frac{k^* b^2 a c}{\mu _A (\gamma + w s)}. \end{aligned}$$
(41)

The biting rate b appears as a square in the definition of \(R_0,\) so b is considered to be an important parameter in determining the prevalence of disease.

From (41), the minimum number of vectors per host necessary to sustain the disease in a control-free environment is found to be

$$\begin{aligned} k_d = \frac{\gamma \mu _A}{b^2 a c} \approx 0.77, \end{aligned}$$
(42)

where default parameter values have been used (Table 1). Hence, even with less than one vector per host, the disease may be maintained at a low level.

B.2 Equilibria and the effect of control

The population and endemic disease equilibria, in the absence of vector control, are found to be

$$\begin{aligned} L^* =&\, \frac{1}{\nu }\left[ \mathrm {e}^{m (\rho / 2 \mu _A - 1) - \mu _L} - 1 \right] ^{\frac{1}{\beta }}, \end{aligned}$$
(43a)
$$\begin{aligned} A^* =&\, \frac{m L^*}{2 \mu _A}, \end{aligned}$$
(43b)
$$\begin{aligned} h^* =&\, \frac{k^* b^2 c a - (\gamma + w s) \mu _A}{k^* b^2 c a + (\gamma + w s) b c}, \end{aligned}$$
(43c)
$$\begin{aligned} v^* =&\, \frac{b c h^*}{b c h^* + \mu _A}. \end{aligned}$$
(43d)

Vector control alters the equilibrium state of the system. Imagine a system in which disease has reached an endemic level [at the level given in (43c)] at some time \(t<t_0,\) and vector control starts at \(t = t_0,\) releasing uA(t) insects (rather than \(u A^*\)) for some constant release ratio u. The new equilibrium at which the system arrives at some time \(t>t_0\) is then of key interest. For early-acting SIT, the new equilibria are

$$\begin{aligned} L^*_E =\, \frac{1}{\nu }\left[ \mathrm {exp}\left\{ \frac{1}{1+u} \frac{m \rho }{2 \mu _A} - m - \mu _L \right\} - 1 \right] ^{\frac{1}{\beta }},\quad A^*_E = \frac{m L^*_E}{2 \mu _A}, \end{aligned}$$
(44)

and (43c) and (43d); and for late-acting SIT the new equilibria are

$$\begin{aligned} L^*_L =\, \frac{1}{\nu }\left[ \mathrm {exp}\left\{ \frac{1}{1+u} \frac{m \rho }{2 \mu _A} - m - \mu _L \right\} - 1 \right] ^{\frac{1}{\beta }},&A^*_L =\, \frac{1}{1+u}\frac{m L^*_E}{2 \mu _A}, \end{aligned}$$
(45a)

and (43c), (43d).

Critical release ratios \(u = u_{\mathrm {c}}\) may be derived which cause disease fadeout (\(h^*\rightarrow 0\)); for early SIT, these release ratios are

$$\begin{aligned} u_{\mathrm {c}} = \frac{\rho m}{2 \mu _A}\left( m + \mu _L + \ln \left[ 1 + \left( 2 (\gamma + w s) \frac{\nu \mu _{A}^{2} N}{b^2 c a m}\right) ^{\beta }\right] \right) ^{-1} - 1, \end{aligned}$$
(46)

while for late SIT \(u_{\mathrm {c}}\) satisfies

$$\begin{aligned} (1 + u_{\mathrm {c}})^{-\beta } \mathrm {exp}\left\{ \frac{1}{1+u_{\mathrm {c}}} \frac{m \rho }{2 \mu _A} - m - \mu _L \right\} = 1 + \left( 2 (\gamma + w s) \frac{\nu \mu _{A}^{2} N}{b^2 c a m}\right) ^{\beta }, \end{aligned}$$
(47)

where solutions of (47) may be found by Newton–Raphson iteration.

Appendix C: A network model

The model outlined above may be extended to allow for movement of mosquitoes between population clusters. Each node (or population cluster) has its own set of state variables \(L_i(t),\) \(A_i(t),\) \(h_i(t),\) \(v_i(t)\) and controls \(u_i(t),\) \(w_i(t)\) satisfying the network model

$$\begin{aligned} \frac{\mathrm {d}L_i}{\mathrm {d}t} =&\, \rho A_i C(A_i, u_i) - f(L_i) L_i - (m + \mu _L) L_i ,\end{aligned}$$
(48a)
$$\begin{aligned} \frac{\mathrm {d}A_i}{\mathrm {d}t} =&\, \frac{m}{2} L_i D(A_i, u_i) - \mu _A A_i - \frac{E_i A_i^{2}}{A_i^*} + \sum _{j\ne i}\sigma _{ij} \frac{E_j A_j^{2}}{A_j^*}, \end{aligned}$$
(48b)
$$\begin{aligned} \frac{\mathrm {d}h_i}{\mathrm {d}t} =&\, \frac{A_i}{N_i} b a (1 - h_i) v_i - (\gamma + w_i s) h_i, \end{aligned}$$
(48c)
$$\begin{aligned} \frac{\mathrm {d}v_i}{\mathrm {d}t} =&\, b c (1 - v_i) h_i - \left( \mu _A + \frac{1}{A_i}\frac{\mathrm {d}A_i}{\mathrm {d}t}\right) v_i - E_i v_i \frac{A_i v_i}{A_i^* v_i^*} + \frac{1}{A_{i}}\sum _{j\ne i}\sigma _{ij} E_j \frac{A_j^{2}v_j^{2}}{A_j^* v_j^*}, \end{aligned}$$
(48d)

where

$$\begin{aligned} C(A_i, u_i) = \left\{ \begin{array}{ll} \displaystyle {\frac{A_i}{A_i + u_i A_i^{*}}, \quad \mathrm {early}\,\, \mathrm {SIT},} \\ \displaystyle {\quad \quad 1, \quad \quad \quad \, \mathrm {late}\,\, \mathrm {SIT},} \end{array} \right.&D(A_i, u_i) = \left\{ \begin{array}{ll} \displaystyle {\quad \quad 1, \quad \quad \quad \, \mathrm {early}\,\, \mathrm {SIT},} \\ \displaystyle {\frac{A_i}{A_i + u_i A_i^{*}}, \quad \mathrm {late}\,\, \mathrm {SIT}.} \end{array} \right. \end{aligned}$$
(49)

The tensor \(\sigma\) defines the weighting of movement between nodes as a function of relative distance. An element \(\sigma _{ij}\) gives the proportion of migrants from the jth population that settle in the ith. For the three-node case, this is defined as

$$\begin{aligned} \sigma = \left( \begin{matrix} 0 &{} \frac{\zeta (d_{12})}{D_2} &{} \frac{\zeta (d_{13})}{D_3} \\ \frac{\zeta (d_{21})}{D_1} &{} 0 &{} \frac{\zeta (d_{23})}{D_3} \\ \frac{\zeta (d_{31})}{D_1} &{} \frac{\zeta (d_{32})}{D_2} &{} 0 \end{matrix}\right) , \end{aligned}$$
(50)

where \(d_{ij}\) is the distance between nodes j and i (thus \(d_{ij} = d_{ji}\)), the function \(\zeta (d_{ij})\) is Laplace-type dispersal kernel,

$$\begin{aligned} \zeta (d_{ij}) = \frac{1}{2}\mathrm {e}^{-|d_{ij}|}, \end{aligned}$$
(51)

and \(D_j = \sum _{k\ne j} \zeta (d_{kj}).\) In this way, \(\sigma\) is normalised such that each column sums to unity, while each row does not. The parameter \(E_i\) is a rate measuring the propensity for emigration of each population cluster. The parameters \(A_i^*\) may be defined as the migration- and control-free equilibria, (43b).

C.1 Optimal control of a network

A new per capita cost functional for a network of n identical population clusters is defined as,

$$\begin{aligned} J_{pc}[\varvec{x}(\varvec{u},\varvec{w})] = \int _0^T \mathrm {e}^{-\psi t} \sum _{i}^{n} \left[ \theta _{1} h_i^{m_{1}} + \theta _{2} h_i w_i^{m_{2}} + \theta _{3} k_i^{*} u_i^{m_{3}}\right] \mathrm {d}t, \end{aligned}$$
(52)

where \(\varvec{u}\) and \(\varvec{w}\) are vectors containing the controls \(u_i\) and \(w_i,\) respectively. The Hamiltonian becomes

$$\begin{aligned} H =&\,\, \mathrm {e}^{-\psi t} \sum _{i}^{n} \left[ \theta _{1} h_i^{m_{1}} + \theta _{2} h_i w_i^{m_{2}} + \theta _{3} k_i^{*} u_i^{m_{3}}\right] + \sum _{i}^{n}\left\{ \right. \lambda _{1}^{(i)}\left[ \rho A_i C(A_i, u_i) - (f(L_i) + m + \mu _L) L_i \right] \nonumber \\&\,+ \lambda _{2}^{(i)}\left[ \frac{m}{2} L_i D(A_i, u_i) - \mu _A A_i - \frac{E_i A_i^{2}}{A_i^*} + \sum _{j\ne i}\sigma _{ij} \frac{E_j A_j^{2}}{A_j^*}\right] \nonumber \\&+ \lambda _{3}^{(i)}\left[ \frac{A_i}{N_i} b a (1 - h_i) v_i - (\gamma + w_i s) h_i\right] \nonumber \\&\,+ \lambda _{4}^{(i)}\left[ b c (1 - v_i) h_i - \left( \mu _A + \frac{1}{A_i}\frac{\mathrm {d}A_i}{\mathrm {d}t} \right) v_i - E_i v_i \frac{A_i v_i}{A_i^* v_i^*} + \frac{1}{A_{i}}\sum _{j\ne i}\sigma _{ij} E_j \frac{A_j^{2}v_j^{2}}{A_j^* v_j^*} \right] \left. \right\} . \end{aligned}$$
(53)

The adjoint system for both the early- and late-acting SIT may be expressed, with the definitions of (49) in mind, as

$$\begin{aligned} \frac{\mathrm {d}\lambda _1^{(i)}}{\mathrm {d}t} =&\, -\frac{\partial H}{\partial L_i} = \lambda _{1}^{(i)}(f'(L_i)L_i + f(L_i) + m + \mu _{L}) - \frac{1}{2} \lambda _{2}^{(i)} m D(A_i, u_i) + \lambda _{4}^{(i)}\frac{m v_i}{2 A_i} D(A_i, u_i), \end{aligned}$$
(54a)
$$\begin{aligned} \frac{\mathrm {d}\lambda _2^{(i)}}{\mathrm {d}t} =&\, -\frac{\partial H}{\partial A_i} = -\lambda _{1}^{(i)}\rho C(A_i, u_i) \left( 1 + A_i \frac{\partial }{\partial A_i}\log {C(A_i, u_i)}\right) - \lambda _{2}^{(i)}\bigg (\frac{m}{2} L_i \frac{\partial }{\partial A_i}D(A_i, u_i) - \mu _{A} \bigg ) \nonumber \\&\, -\lambda _{3}^{(i)}\frac{b a}{N_i}(1 - h_i) v_i - \lambda _{4}^{(i)} v_i \Bigg (\frac{m L_i}{2 A_i^{2}} D(A_i, u_i) \left[ 1 - A_i \frac{\partial }{\partial A_i}\log {D(A_i, u_i)} \right] - \frac{E_i v_i}{A_i^* v_i^*}\Bigg ) \nonumber \\&\, + \lambda _{4}^{(i)}\frac{1}{A_i^2}\sum _{j \ne i} \sigma _{ij} E_j \frac{A_j^{2}v_j^{2}}{A_j^* v_j^*} - 2 \frac{E_i A_i v_i^2}{A_i^* v_i^*} \sum _{j \ne i} \frac{\lambda _4^{(j)}}{A_j}\sigma _{ji} + 2 \frac{E_i A_i}{A_i^*}\left( \lambda _2^{(i)} - \sum _{j \ne i} \lambda _2^{(j)} \sigma _{ji} \right) , \end{aligned}$$
(54b)
$$\begin{aligned} \frac{\mathrm {d}\lambda _3^{(i)}}{\mathrm {d}t} =&\, -\frac{\partial H}{\partial h_i} = -\mathrm {e}^{-\psi t}(m_{1} \theta _{1} h_i^{m_{1}-1} + \theta _{2} w_i^{m_{2}}) + \lambda _{3}^{(i)}\left( \frac{A_i}{N_i} b a v_i + \gamma + w_i s\right) - \lambda _{4}^{(i)} b c (1 - v_i), \end{aligned}$$
(54c)
$$\begin{aligned} \frac{\mathrm {d}\lambda _4^{(i)}}{\mathrm {d}t} =&\, -\frac{\partial H}{\partial v_i} = -\lambda _{3}^{(i)}\frac{A_i}{N_i} b a (1 - h_i) + \lambda _{4}^{(i)}\left( b c h_i + \mu _A + \frac{1}{A_i}\frac{\mathrm {d}A_i}{\mathrm {d}t}\right) \nonumber \\&+ 2 E_i \frac{A_i^2 v_i}{A_i^* v_i^*}\left( \frac{\lambda _4^{(i)}}{A_i} - \sum _{j \ne i} \frac{\lambda _4^{(j)}}{A_j} \sigma _{ji}\right) , \end{aligned}$$
(54d)

where \(1 - \sum _{j\ne i}\sigma _{ji} = 0\) has been used in (54b).

The optimality conditions that define the best-practice for using artemisinin and vector-release controls are

$$\begin{aligned} \frac{\partial H}{\partial u_i} = 0 \quad \mathrm {at} \,\, u_i = u_i^{*},\quad \frac{\partial H}{\partial w_i} = 0 \quad \mathrm {at} \,\, w_i = w_i^{*}. \end{aligned}$$
(55)

The controls \(u_i\) and \(w_i\) never appear in equations \(\mathcal {L}_j \mathbf{x}_j,\) for nodes i, j, differential operator \(\mathcal {L}\) and vector of states and adjoints \(\mathbf{x}\)—in either (48) or (54)—where \(j \ne i,\) hence the results of the previous section are still valid for the network model. Then, the optimal controls for each pricing structure and each population cluster may be defined as in (18), (21) and (29).

Rights and permissions

Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Khamis, D., El Mouden, C., Kura, K. et al. Optimal control of malaria: combining vector interventions and drug therapies. Malar J 17, 174 (2018). https://0-doi-org.brum.beds.ac.uk/10.1186/s12936-018-2321-6

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://0-doi-org.brum.beds.ac.uk/10.1186/s12936-018-2321-6

Keywords