Royal Society Publishing

Slowing the evolution of insecticide resistance in mosquitoes: a mathematical model

Stephen A. Gourley, Rongsong Liu, Jianhong Wu


A big problem in malaria control is the rapidity with which mosquitoes can develop resistance to insecticides. The possibility of creating evolution-proof insecticides is therefore of considerable interest. Biologists have suggested that effective malaria control, with only weak selection for insecticide resistance, could be achieved if insecticides target only old mosquitoes that have already laid most of their eggs. The strategy aims to exploit the fact that most malarial mosquitoes do not live long enough to transmit the disease. We derive, analyse and compare two mathematical models, one for an insecticide that kills on exposure, and the other for an insecticide that targets only older mosquitoes. Both models predict that insecticide-resistant mosquitoes will become dominant over time but, very importantly, this occurs on a very much slower time scale when the insecticide only affects older mosquitoes. We present analytical results on linear and global stability of the non-trivial equilibrium in which only the resistant mosquito strain is present, together with a theorem comparing the rates of convergence for the two models. Numerical simulations show that the effect of targeting only old mosquitoes on the evolution of resistance is dramatic.

1. Introduction

The idea of evolution-proof insecticides is currently a hot topic in mosquito control (Read et al. 2009) because of the potential to slow down or even halt the evolution of insecticide resistance in mosquitoes. While insecticides are an effective and cheap method of controlling malaria, it is well known that mosquitoes can develop resistance and that this can happen on a surprisingly fast time scale. It was known by the end of the 1960s that the mosquito species that transmit malaria were developing resistance to dichlorodiphenyltrichloroethane (DDT), the single insecticide that had been relied on until then (see Hemingway et al. (2002) and Kelly-Hope et al. (2008) for recent discussions on the management of insecticide resistance). Current insecticides kill extremely rapidly after contact, but their high lethality leads to intense selection for resistance because they kill young female adults. The very limited insecticide arsenal available (there are still just four classes, only one of which is approved for use on bed nets) is leading to an increasing focus on resistant-management strategies, whereby the useful lifespan of an existing insecticide is enhanced by the way in which the insecticide is used. Resistant-management strategies may include using different insecticides in space and time and restricting use to specific times and locations. Indeed, it seems clear that the development of new insecticides will not in itself provide a sustainable solution to the problem of resistance. All existing insecticides were new at some point, development costs are very high and mosquito evolution dooms them all to failure. See Koella et al. (2009) where it is further suggested that, since the Second World War, the operational life of an insecticide, in areas of very widespread use, has been of the order of just 5 years.

Very importantly, it is known that most mosquitoes do not actually live long enough to be able to transmit malaria. This is due to a relatively long latency stage during which the malaria parasite has to go through various developmental stages and replicative cycles in the mosquito before reaching the salivary glands from which the sporozoites can be transmitted via bites to humans. The duration of the latency period is of the order of 10–14 days in areas of high malaria transmission (Charlwood et al. 1997; Killeen et al. 2000). This is almost a lifetime to a mosquito, and therefore the majority of eggs produced by a female mosquito in her lifetime are laid before the mosquito is able to transmit malaria to a human. To exploit these facts to advantage, Read et al. (2009) propose that the useful lifespan of an insecticide can be enhanced if the insecticide targets only old mosquitoes. They call such an insecticide a late-life acting (LLA) insecticide. More precisely, the insecticide should start to take effect on a female mosquito after she has laid most of her eggs but before the mosquito can start to transmit malaria. This results in much weaker selection for resistance. Koella et al. (2009) suggest that, in areas where DDT control worked well for 5 years with high and continual coverage, an insecticide with a delayed action could provide continuous control for as long as 35 years.

There are two ways in which an insecticide might kill only old mosquitoes and not younger adult mosquitoes. One possibility is that the insecticide could begin to take effect not on initial contact but some time later. This could perhaps be achieved by cumulative exposure over time to doses each of which would not be lethal. Existing insecticides could be used with doses lower than those currently in use (or the insecticide could simply be diluted). An alternative is the use of fungal pesticides that kill a mosquito 7–14 days after contact (see Thomas & Read 2007 and ongoing work by these investigators). The second approach is to have an insecticide that takes effect only on older insects by exploiting physiological changes and weaknesses associated with ageing.

In this paper, we consider the first of the two possibilities mentioned above, i.e. the insecticide does not take effect on first contact but only after a time delay. We propose to divide the lifespan of the female mosquito into three stages: the larval stage (which is understood to include all the pre-adult stages of egg, larva and pupa and is of total duration τi), the adult stage and the ‘old age’ stage. We assume that once a mosquito matures from a larva into an adult, it is immediately exposed to insecticide. Such exposure may not in reality be immediate, but an adult female needs to find a blood meal before it can lay eggs and the search for such a blood meal is likely, in areas of intensive mosquito control using insecticides, to bring the female mosquitoes into contact with insecticides on bed nets or on the walls of homes. As an approximation it is assumed that those mosquitoes contact the insecticide immediately on maturation and then there is a delay, of duration τa, before the insecticide takes effect. Once this starts to happen the mosquito is classed as an old mosquito. Thus ‘old age’, in this paper, does not necessarily mean that the mosquito is near the end of its life or has undergone particular physiological changes owing to ageing. It simply means that the mosquito has been an adult for at least τa units of time and the insecticide has started to take effect. The parameter τa is within our control, if we can control the length of the delay before the insecticide takes effect.

We consider two mathematical models. The first is for the case of an insecticide that kills on exposure, and not after a delay, and we refer to it as the conventional insecticide model. The adult mosquitoes are classed as V (vulnerable) or R (resistant) according to whether the insecticide can have an effect or not. Even this simple model involves time delay owing to the fact that, at a given time t, the rate at which larvae mature into adults is related to the birth rate at the earlier time tτi, and therefore to the number of egg-laying adults at that time. For this model we establish, in §2a, the local and global stability of the equilibrium in which the resistant strain is established and the vulnerable strain is extinct. Then, in §2b, we develop a more complex model for the case of an insecticide that targets only older mosquitoes of age exceeding τi+τa. For this model, the linear stability analysis about the non-trivial equilibrium in which the vulnerables are extinct is tricky but tractable. A comparison of the dominant eigenvalues for the two models is carried out and we prove that, in the model for an insecticide with delayed action, the resistant strain still becomes the dominant species but on a slower time scale than in the case when the insecticide acts instantly. As is usual in time-delay systems, the characteristic equations are transcendental, but numerical computations of the dominant eigenvalues demonstrate that the prolongation of the evolution of resistant mosquitoes can be dramatic if the insecticide acts with a delay. Neither of the models takes account of fitness costs of resistance. Inclusion of costs could alter the outcomes and this point will be discussed later.

2. Derivation and analysis of the models

(a) Insecticide that kills on exposure

We let V (t) denote the number of female mosquitoes that are vulnerable to the insecticide and R(t) the number that are resistant to it, and we derive a system of delay differential equations to be solved for V (t) and R(t). Resistance arises via a genetic mutation that generates a small non-zero initial value for R(t). The model derivation proceeds by reducing an age-structured model, for the densities v(t,a) and r(t,a) of vulnerable and resistant mosquitoes at time t of age a, to a system of delay equations. It is similar to the derivation of the more complex model to be considered in §2b, so here we only present brief details. If insecticide does not affect larval mosquitoes then we may assume that the density v(t,a) of vulnerables satisfies Embedded Image2.1where the density v(t,a) is defined such that for infinitesimal da the number of vulnerables aged between a and a+da is v(t,a) da. In equation (2.1), μi is the per capita natural death rate for immature (larval) mosquitoes and τi is the maturation time for the larva from egg to adult and is the total duration of all pre-adult stages. For adult mosquitoes, we let their natural per capita death rate be μa, and we let μa+δ denote the overall per capita death rate for vulnerables (so that δ is the per capita insecticide-induced mortality). Therefore, for adult vulnerables, Embedded Image2.2The total number of adult vulnerables is Embedded Imageand, using equation (2.2) and assuming that Embedded Image, it satisfies Embedded Image2.3We calculate v(t,τi) from equation (2.1) by integrating along the characteristics, obtaining v(t,τi)=v(tτi,0)eμiτi. In this type of model formulation, v(t,a) has an additional interpretation: it is the rate at time t at which individuals pass through age a, and therefore v(t,0) is the birth rate (of vulnerables). We assume that vulnerable and resistant mosquitoes of all adult ages are equally likely to lay eggs, that egg production is continuous and that all adult mosquitoes compete with each other on an equal basis for resources. Thus, the overall egg-laying rate at time t will be taken to be a function b(M(t)) of the total number of adult mosquitoes, M(t)=V (t)+R(t). On the assumption that insecticide resistance is an inherited characteristic involving changes in insect genes (Hemingway et al. 2004), we assume that the offspring of vulnerables are vulnerable and the offspring of resistant mosquitoes are resistant, so that the birth rates for vulnerable and resistant mosquitoes are given, respectively, by Embedded Image2.4Thus, Embedded Image, and so equation (2.3) becomes the V equation of the system Embedded Image2.5The R equation of system (2.5), for resistant mosquitoes, is derived in very much the same way, except that there is no δ term as these mosquitoes are resistant to insecticide.

We always assume that b(0)=0 and that Embedded Image. System (2.5) therefore has no singularity and always has (0,0) as an equilibrium. We show under some reasonable assumptions on the function b(⋅) that, according to model (2.5), if δ>0, the vulnerable mosquitoes are driven to extinction while the resistant ones take over. What is particularly of interest, in both these sections and the model of §2b, is the time scale on which this process takes place. This can be estimated from linearized analysis at the relevant equilibrium. It is straightforward to see that system (2.5) has no equilibrium with both V,R>0. Apart from the zero equilibrium, it may have an equilibrium with V =0 and R>0 and another with V >0 and R=0. We prove that the former is locally stable in theorem 2.1. In this equilibrium, V =0 and R=R*, where R* satisfies Embedded Image2.6If an equilibrium (V,R)=(0,R*) exists, then another equilibrium (V*,0) also exists if δ is sufficiently small, and analysis similar to the proof of theorem 2.1 shows it to be linearly unstable. Such an equilibrium does not exist if δ is large (i.e. the insecticide is very effective). Small δ is arguably more realistic, as a single mutation is likely to confer only a small increase in insecticide tolerance.

Theorem 2.1

Suppose that δ>0, that R*>0 satisfies equation (2.6) and that μa>eμiτib′(R*)>0. Then, the equilibrium (V,R)=(0,R*) of system (2.5) is locally asymptotically stable.


Linearization of system (2.5) at the equilibrium (0,R*) yields the partially decoupled system Embedded Image2.7where Embedded Image. In view of equation (2.6), the first of these equations can be put in the form Embedded Image2.8and, since δ>0, Kuang (1993, theorem 3.2.1) yields that Embedded Image as Embedded Image. The Embedded Image equation then becomes Embedded ImageSince μa>eμiτib′(R*)>0, Embedded Image as Embedded Image (Kuang 1993, theorem 3.2.1). ■

For the global stability of (0,R*), to be proved next, we shall need to assume that the birth function b(⋅) has the classic hump shape with the equilibrium on the increasing side, in the following precise sense: Embedded Image2.9The Nicholson’s blowflies birthrate, b(M)=pMeqM, satisfies equation (2.9) for appropriate parameter values.

Theorem 2.2

Suppose that δ>0, and R* satisfies equation (2.6) and that b(⋅) satisfies assumption (2.9). Then, the equilibrium (V,R)=(0,R*) of system (2.5) is globally attractive for all solutions with non-negative initial data, such that R(θ)≢0 on [−τi,0].


First note that V (t)≥0 and R(t)≥0 for all t>0; this follows from Smith (1995, theorem 5.2.1 on p. 81). Next, let us prove that both V (t) and R(t) enter, in finite time, the interval of values where b(⋅) is increasing. To see this, note that, if we add the equations in system (2.5), Embedded Image2.10so that Embedded Image. A graphical argument using the properties of b(⋅) shows that eμiτibmax/μa<Rmax. Therefore, M(t)<Rmax for sufficiently large t, and the same is true for V (t) and R(t). Inequality (2.10) can also be used to prove the positive invariance of the region Embedded ImageWe proceed via the use of theorem B in Hsu et al. (1996). It is a rather general result that admits infinite-dimensional dynamical systems that have solutions which evolve in a product of two partially ordered Banach spaces. Each space contains the solutions involving one component (V or R, in our case) remaining identically zero. If the solution operator is strictly order preserving, the origin is a repelling equilibrium, each subsystem with V (t)≡0 or R(t)≡0 contains a global attractor, and some other technical conditions are satisfied, then, if there is no coexistence equilibrium, solutions must approach one of the boundary equilibria. Linear stability theory can be used to decide which boundary equilibrium is approached. The theorem is stated in the entire phase space C([−τi,0],R2), but it remains valid when applied to the region Embedded Image owing to the attractivity and invariance of the region within which b(⋅) is increasing. So, in what follows, we assume b′>0 without loss of generality.

How the proof proceeds depends on whether or not there is an equilibrium with V =V* and R=0. For such an equilibrium, V* would satisfy (μa+δ)V*=eμiτib(V*). Under hypothesis (2.9), such a boundary equilibrium exists if δ is sufficiently small and we assume this to be so. However, the theorem remains true, even if δ is large enough to preclude the existence of such an equilibrium.

Let X=C([−τi,0],R) and let X+ be the positive cone of X. Define the cone K=X+×(−X+), which induces the partial ordering Embedded Imagewhere ≤ is the partial ordering induced by X+. To verify hypothesis (H1) in Hsu et al. (1996, p. 4086), it is necessary to check that the solution operator T(t), which maps the initial state (V (θ),R(θ)), θ∈[−τi,0] to the state at time t, namely (Vt(θ),Rt(θ)), θ∈[−τi,0], is strictly order preserving with respect to <K. Here, < means ≤ and ≠, and Vt(⋅) is the function Vt(θ)=V (t+θ), θ∈[−τi,0]. So we need to check that Embedded Image implies that Embedded Image, i.e. that Embedded Imagewith Embedded Image or Embedded Image on [−τi,0], implies Embedded Imagewith Embedded Image or Embedded Image, θ∈[−τi,0]. This amounts to showing strict order preserving for the variable (V,−R), in the sense of the usual partial ordering in R2. Letting Embedded Image and Embedded Image, we find the evolution equations for Embedded Image to be Embedded Image2.11where Embedded ImageRoutine differentiation, using that b is increasing and that Embedded Image, shows that F1 and F2 are each increasing with respect to both variables if b′(M)<b(M)/M, and this is easily seen to follow from the hypothesis that b′′<0 on the increasing side of b. Therefore, theorem 5.1.1 (p. 78) of Smith (1995) yields that system (2.11) preserves ordering. To show the order preserving is strict, suppose that, for some t*>0, Embedded Image and Embedded Image for θ∈[−τi,0]. Then, the two solutions (V,R) and Embedded Image of system (2.5) agree on an interval of length τi and, by uniqueness, Embedded Image for all t>t*. We show that this contradicts Embedded Image or Embedded Image on [−τi,0] by showing that it implies Embedded Image and Embedded Image on [−τi,0]. Introducing Embedded Image and Embedded Image, these functions are zero for ti, where N is an integer such that it*. Hence, from system (2.5), Embedded Imagefor t≥(N−1)τi, where Embedded ImageUnder the hypotheses, G1 is increasing with respect to V and decreasing with respect to R, while G2 is increasing with respect to R and decreasing with respect to V . Therefore, for t≥(N−1)τi, Embedded Imagefor some function Embedded Image. Since ∂G1/∂V >0, this implies ξ(t)≤0, and so ξ(t)=0. We have shown that ξ(t)=0 for ti implies the same for t≥(N−1)τi. The argument can be continued to yield ξ(t)=0 for t∈[−τi,0], and the same can be shown for η(t) using properties of G2. So, Embedded Image and Embedded Image on [−τi,0], giving a contradiction. This establishes (H1) in Hsu et al. (1996).

To verify (H2) in Hsu et al. (1996), we need to check that (0,0) is a repelling equilibrium. This is easy to do and follows from μa<eμiτib′(0). Hypothesis (H3) in Hsu et al. (1996) holds because {0}×X+ is invariant (if V starts zero, it stays zero) and since the equilibrium R* of R′(t)=−μaR(t)+eμiτib(R(tτi)) is globally attracting within {0}×X+ under our hypotheses (see Kuang 1993, theorem 4.9.4, p. 164). Also, the equilibrium (V,R)=(V*,0), if δ is small enough to ensure its existence, is globally attracting within X+×{0}.

To verify (H4) of Hsu et al. (1996), note that if both V and R start off not identically zero, they both become and remain strictly positive. Verification of (H4) also involves checking that if Embedded Image and either (V (⋅),R(⋅)) or Embedded Image belongs to int (X+×X+), then Embedded Image for t>0. In terms of the functions ξ(t) and η(t), we need to show that if ξ(θ)≥0 and η(θ)≥0 on [−τi,0], ξ(θ)≢0 or η(θ)≢0, then ξ(t) and η(t) become and remain strictly positive. The evolution equation for ξ(t), in the form Embedded Imageyields Embedded Image

The requirement that (V (⋅),R(⋅)) or Embedded Image belongs to int (X+×X+) guarantees that φ(tτi) and Embedded Image cannot both be zero, so that Embedded Image. Hence, if ξ(θ)≢0 on [−τi,0], then the usual argument (of considering the initial evolution in the interval t∈[0,τi]) yields that ξ(t) must become positive at some time in [0,τi]. Having become positive, ξ(t) remains positive because ξ′(t)≥−(μa+δ)ξ(t). The possibility that η(θ)≢0 on [−τi,0] is treated similarly by making use of the evolution equation for η(t).

Having checked these hypotheses, the trichotomy in Hsu et al. (1996, theorem B) applies and, in the absence of a coexistence equilibrium, we conclude that solutions must approach one of the boundary equilibria (V*,0) or (0,R*). The former can easily be shown to be linearly unstable to perturbations involving the introduction of R (the linearized equation is R′(t)=−μaR(t)+(μa+δ)R(tτi)) and the latter is linearly stable by theorem 2.1. The proof is complete. ■

Later in the paper, we will be particularly interested in the decay rate λall, which tells us something about how fast the solutions evolve to the equilibrium (0,R*), at which vulnerables are extinct, for the present model in which all adult mosquitoes, and not just older ones, are targeted. The parameter λall is the dominant root (the root of the largest real part) of the characteristic equation Embedded Imagewhich results from the use of the trial solution Embedded Image in equation (2.8). The dominant root is real and negative (Smith 1995, theorem 5.5.1, p. 92) and the time taken to reach equilibrium is of the order |λall|−1. Later, we will also be concerned with the corresponding decay rate λold associated with the model of the next section, in which only older mosquitoes are targeted.

(b) Insecticide that targets older mosquitoes

Again V and R denote, respectively, the numbers of female adult vulnerable and insecticide-resistant mosquitoes, but in this section mosquitoes are considered to have three stages of life in all. These are (i) the larval stage, considered to be of duration τi, (ii) the adult stage, which is of duration τa and does not include old age, and (iii) the ‘old age’ stage, for mosquitoes of age exceeding τi+τa. Subscripts ‘i’, ‘a’ and ‘o’, standing for immature, adult and old, are used as appropriate so that, for example, Va is the number of adult vulnerable mosquitoes and Vo is the number of old vulnerable mosquitoes. We propose the following model, which we call the LLA insecticide model: Embedded Image2.12where Ma(t)=Va(t)+Ra(t), etc. Note that in this system, δ, which represents insecticide-induced death, only appears in the equation for old vulnerable mosquitoes, and not all vulnerable mosquitoes. This is how we model the age-dependent effect of the insecticide. Old mosquitoes may also have a different natural per capita mortality μo from the corresponding mortality μa for adults. In reality, one expects that μoμa. Let us derive in detail the Va equation; the derivations of the other equations are similar. Adults are individuals aged between τi and τi+τa and therefore the number of vulnerable adults is Embedded Image2.13For vulnerable mosquitoes at all stages of life, their age density v(t,a) satisfies Embedded Image2.14where Embedded Image2.15Differentiating equation (2.13) and using equations (2.14) and (2.15) gives Embedded Image2.16The last two terms of this are the maturation rate and the rate of reaching old age, and these two terms can be calculated in terms of the birth rate v(t,0). Solving equation (2.14) for a general μ(a) yields that Embedded Image

With μ(a) given by expression (2.15), this gives Embedded Image2.17and Embedded Image2.18In the second expression, the exponential factors represent the probabilities of surviving the larval phase and the adult phase to reach old age. The birth rate for vulnerables is v(t,0). Adult and old, vulnerable and resistant mosquitoes are assumed to be equally likely to lay eggs, so that the overall egg-laying rate is b(Ma(t)+Mo(t)). As in the previous section, we deduce that Embedded Image2.19because the above expression represents the proportion of the eggs laid that contain the genes for vulnerability, having acquired those genes from their parents. From equation (2.19), we calculate (2.17) and (2.18) and insert the results into equation (2.16), completing the derivation of the Va equation of system (2.12).

The Va equation of system (2.12) essentially formulates the change rate of the Va population as the combination of death, the maturation rate, which is the birth rate at time tτi, and the progression into the ‘old age’ class, which happens τa time units after maturation. Integration of this equation yields the integral equation Embedded Imageif this matching condition is met initially. So system (2.12) can be replaced by the differential equations for Vo and Ro, coupled with the above integral equation and its analogue for Ra. In what follows, we shall not distinguish the system (2.12) from its equivalent formulation involving integral equations as just described.

As in the conventional insecticide model of §2a, we shall be concerned primarily with the linear stability of the equilibrium in which the vulnerables are extinct, but the resistant mosquitoes are not. So, at equilibrium, Va=Vo=0 and Embedded Image, Ro=R*o where Embedded Image2.20Thus, Embedded Image and we have a single equation for R*o, Embedded Image2.21

This makes it easy to impose biologically realistic conditions sufficient for the existence of unique equilibrium values with R*o>0 and R*a>0. Such a reasonable assumption would be Embedded Image2.22For the linear stability of the equilibrium Embedded Image of system (2.12), we shall need the following additional condition: Embedded Image2.23We now investigate this linear stability. Near the equilibrium Embedded Image, the linearization of the V equations is Embedded Image2.24Setting Embedded Image and Embedded Image yields the linearization of the R equations, Embedded Image2.25and Embedded Image2.26System (2.24) involves only the V variables and therefore the linearized system is partially decoupled. The characteristic equation, which we will not write down in full, decomposes into a product of factors leading effectively to two characteristic equations, both of which must be considered. The first of these relates to system (2.24), and to enable satisfactory analytical progress, it is essential that this characteristic equation be cast into a suitable form. Note that the first equation of system (2.24) can be recast as Embedded Image2.27We solve equation (2.27) together with the second of system (2.24) using the usual ansatz Embedded Image, where c1 and c2 are constants. After some algebra, this leads to the following useful form for the characteristic equation to be solved for λ: Embedded Image2.28where Embedded Image2.29In theorem 2.3, we prove that all roots of equation (2.28) have a negative real part, and this means that Embedded Image as solutions of system (2.24). But this means that in the study of the linearization of the R equations, it is sufficient to set the V variables to zero, leading to the simpler linearized equations Embedded Image2.30the first of which can be rewritten as Embedded Image2.31This leads to the second characteristic equation Embedded Image2.32where Embedded Image2.33

We prove the following result on the linear stability of Embedded Image.

Theorem 2.3

Suppose that assumption (2.22) and inequality (2.23) hold. Then, for any δ>0, the equilibrium Embedded Image of system (2.12) is locally asymptotically stable.


We prove that all roots of equation (2.28) satisfy Re λ<0. The situation for equation (2.32) is similar. Let us first prove that Embedded Image2.34Writing Embedded Image, it is necessary to show that, for x≥0, Embedded Image2.35or, equivalently, that Embedded Image2.36where Embedded Image. Let us first prove that, for a fixed Embedded Image, Embedded Image2.37Since Embedded Image, it is clearly enough to check that h(ξ,η*)<h(ξ,0) where, for a fixed ξ, η* is any turning point of the function h(ξ,⋅) other than zero. Such a turning point satisfies Embedded Imageso that Embedded ImageTo show h(ξ,η*)<h(ξ,0), it is sufficient to check that eξ<(1−eξ)2/ξ2. This is true for all ξ>0; indeed, it is equivalent to the assertion that j(ξ)>0, where j(ξ)=1−eξξeξ/2. But this function satisfies j(0)=0 and Embedded Image for ξ>0, so that j(ξ)>0 when ξ>0. Having established equation (2.37), to show (2.36) it is necessary to show that Embedded Image2.38and this follows from the fact that the function Embedded Image is decreasing (see Gourley et al. 2007). We have shown inequality (2.34).

Now suppose there exists a root Embedded Image of equation (2.28) such that Embedded Image. Then, Embedded Image2.39Writing f(λ) in the form Embedded Image2.40where g(x)=(1−ex)/x, we have, since Embedded Image, Embedded ImageIt is our claim that, for Embedded Image, Embedded Image2.41The right-hand side of this simplifies to Embedded Image because of (2.20). Therefore, in view of the inequality |z1z2|≥|z1|−|z2| for z1,z2C, it suffices to show that Embedded Image2.42But, using the first equation of (2.20) and the definition of g, we can show that inequality (2.34) is equivalent to Embedded Imageand inequality (2.42) follows immediately, since Embedded Image. Using inequality (2.41) and the simplification of its right-hand side, we may now deduce from equation (2.39) that Embedded Imageusing the second equation of (2.20). But for Embedded Image, this is impossible because it implies that Embedded Image is within a disc contained entirely in the open left-half plane.

We now treat equation (2.32). The above argument essentially proves that, if conditions are such that the dominant real root of equation (2.28) is negative, then there can be no complex roots with positive real part. The arguments for equation (2.32) are similar, so for simplicity we only look at its real roots. Simple graphical arguments yield that the dominant real root of equation (2.32) is negative if we can show that Embedded Image2.43and f1(λ) is decreasing for λ≥0. Note that inequality (2.23) implies Embedded Imageso that Embedded Image2.44Since inequalities (2.23) and (2.44) hold, we have Embedded Imageso inequality (2.43) holds. To show that f1(λ) is decreasing for λ≥0, we can write it in terms of the positive and decreasing function g, similar to expression (2.40) for f(λ). Written in this form, f1(λ) has a decreasing numerator and increasing denominator (note that Embedded Image, from inequality (2.23)). It follows that f1(λ) is decreasing if Embedded Imagewhen λ≥0. The left-hand side is increasing in λ and so it is enough to check the inequality for λ=0, but in that case, it reduces to inequality (2.44). ■

In the next theorem, we prove that, if μo=μa, then, although the solutions of models (2.5) and (2.12) both converge to equilibria in which the resistant mosquitoes are present and the vulnerables are extinct, in the latter model, the convergence to equilibrium occurs on a slower time scale. Insecticide resistance therefore develops more slowly if the insecticide affects only the older mosquitoes. Later, we present numerical computations which demonstrate that the effect is dramatic.

Theorem 2.4

Suppose that μoa and δ>0. Then, the convergence of solutions to the equilibrium Embedded Image of system (2.12) occurs more slowly than the convergence of solutions of system (2.5) to its equilibrium (0,R*). More precisely, if λold and λall are, respectively, the dominant real eigenvalues of the linearizations of systems (2.12) and (2.5) at these equilibria, then Embedded Image2.45


The eigenvalue λall satisfies Embedded Image2.46and, since μo=μa, λold satisfies equation (2.28), i.e. Embedded Image2.47Define f1(λ)=μaeλτi and Embedded Imagewith f(λ) given by (2.29). Sketches of the curves Embedded Image, Embedded Image and Embedded Image and continuity arguments show that, if Embedded Image2.48then inequality (2.45) follows. Since we assume μo=μa, we may add the equations in (2.20) to get Embedded Image, and so Embedded ImageAlso, f1(λall)=μaeλallτi. So, to show that f2(λall)>f1(λall), it is necessary to show that Embedded Imagewhich is equivalent to Embedded Imageand the truth of this follows from equation (2.46). It is easily checked that f2(0)=μa<μa+δ. Thus, (2.48) holds and the proof is complete. ■

3. Numerical simulations and discussion

In this section, we use numerical simulation to confirm the analytical results and explore other properties of the models. Figures 1 and 2 illustrate the two kinds of behaviours of models (2.5) and (2.12). We chose the egg-laying rate as b(M)=pMeqM, which is the famous Nicholson’s blowflies birthrate (Gurney et al. 1980). It is a popular choice in modelling insect dynamics owing to its positivity, the ease of interpretation of the parameter p (per capita egg production rate at low densities) and the fact that it models reduction in egg production at very high densities owing to intraspecific competition for resources. If p=15, both systems evolve to a state in which the vulnerable mosquitoes have died out owing to the insecticide, and the resistant mosquitoes have become the dominant species with their numbers stabilized. Figure 1a shows the evolution for the conventional insecticide model (2.5) in which the insecticide kills mosquitoes on exposure, and not after a delay. This is qualitatively very similar to Koella et al. (2009, fig. 1b), in fact, the time scales on which the resistant individuals take over is similar. The main notable difference between the evolution of solutions of the two models (2.5) and (2.12) is, of course, the time scale on which the resistant strain takes over. It is very much slower in the case of an insecticide that only targets old mosquitoes. Note, in the latter case, that Ra+Ro actually evolves to a lower equilibrium value than does R in system (2.5). This is because, in this simulation, we use parameter values with μo>μa since older mosquitoes have a higher per capita natural mortality.

Figure 1.

In this figure, p=15. (a) The evolution for the conventional insecticide model (2.5). For this model, the vulnerable mosquito dies out and the resistant strain very rapidly takes over, the process being essentially complete after around 250–300 days. Dot-dashed line, V (t); solid line, R(t). The culling rate δ=0.4. For comparison, (b) shows the evolution for the LLA insecticide model (2.12), again with δ=0.4, but this time the insecticide only affects older vulnerable mosquitoes. Here, the vulnerables are again driven to extinction by the resistant strain, but on a much longer time scale of the order of 700 days. Dot-dashed line, Va(t)+Vo(t); solid line, Ra(t)+Ro(t). Other parameter values are taken from table 1.

Figure 2.

In this figure, p has been increased to p=50 and the death rate of larvae μi decreased to μi=0.25, with other parameter values being the same as in figure 1. (a) The evolution for the conventional insecticide model (2.5): the vulnerable mosquitoes have died out and the resistant strain has almost completely taken over by 150 days, with the culling rate δ=0.4. Dot-dashed line, V (t); solid line, R(t). (b) Still with δ=0.4, this figure shows the corresponding evolution for the LLA insecticide model (2.12). Dot-dashed line, Va(t)+Vo(t); solid line, Ra(t)+Ro(t). The effect is a substantial lengthening, to around 600 days, of the time taken for the resistant strain to take over. In contrast to figure 1, the numbers of resistant mosquitoes evolve to a limit cycle rather than an equilibrium. This change is attributable to the increase in the value of p and decrease in the value of μi. Since LLA insecticides are not yet in use, it is not known whether natural mosquito populations can oscillate as a result of being subjected to an LLA insecticide.

View this table:
Table 1.

Parameter values used, except where figure captions state otherwise. The values relate mainly to the malarial Anopheles genus.

Figure 2 shows that, if we increase p to 50, even though the vulnerable mosquitoes still die out owing to the insecticide and the resistant mosquitoes still become dominant, in model (2.12), the resistant mosquitoes evolve to an oscillatory state rather than an equilibrium. Note that, when p is large, the R value of the equilibrium being considered is beyond the interval in which the birth rate function b(⋅) is monotone increasing.

Figures 3 and 4 demonstrate how the time to extinction of vulnerables (defined as the reciprocal of the modulus of the dominant eigenvalue) for each model changes with respect to two of the parameters, which we consider to be particularly important owing to the ability to change them as part of a control strategy: the insecticide-induced death rate, δ (which is effectively a measure of the effectiveness of the insecticide) and the duration τa of the adult stage of the mosquito. It is very important to stress here that the duration of the adult stage is defined solely by the action of the insecticide, and not by anything intrinsic to the mosquito such as the onset of particular physiological changes attributable to age. It is assumed that mosquitoes become exposed to the insecticide on maturation, but that the action of the insecticide is delayed. What we call ‘old age’ is just the age, namely τi+τa, at which the insecticide starts to take effect on the mosquito. A mosquito is classed as an adult if it has matured, and has been exposed to the insecticide since then, but the insecticide has not yet started to take effect on the mosquito. In this sense, the duration τa of the adult phase is within our control if we can design the insecticide to act with a specified time delay τa.

Figure 3.

In this figure, μo=0.1, the same as μa and p=50. (a,b) Show how the time to extinction of vulnerables (defined as 1/|λall| and 1/|λold| for models (2.5) and (2.12), respectively) changes as we increase the culling rate δ or the duration τa of the adult stage. The plots confirm inequality (2.45) and also yield quantitative information. Increasing δ (essentially, the effectiveness of the insecticide) shortens the time taken for the resistant strain to take over. Very importantly, this effect is many times more dramatic in the conventional insecticide model (2.5) than in the case where only older mosquitoes are targeted. (b) Shows that if we increase the delay τa before the insecticide takes effect, then the time for the resistant strain to take over increases, and the effect is significant. The conventional insecticide model does not have a parameter τa, hence the horizontal line. (a) τa=10. (b) δ=0.4. (a,b) Dot-dashed line, conventional insecticide model; solid line, LLA insecticide model.

Figure 4.

(a,b) In this figure, μo=0.25 and other parameters are the same as in figure 3. We draw similar conclusions to those described in the caption to figure 3, but here we illustrate a situation when μaμo. Dot-dashed line, conventional insecticide model; solid line, LLA insecticide model.

From figures 3 and 4, we see that, whether μa equals μo or not, for the conventional insecticide model (2.5), vulnerables go extinct more quickly than in the LLA insecticide model (2.12). See the caption of figure 3 for further comments.

 Read et al. (2009) point out that the evolution of resistance to an LLA insecticide could be slowed down even more if the insecticide could be made to act only on malaria-infected mosquitoes because, in this way, one further relaxes selection for resistance without any loss of control. Moreover, there would be increased selection pressure favouring mosquitoes that are resistant to malaria.

Another aspect that we have not emphasized in this paper is the possibility of larviciding (killing larval mosquitoes using larvicides). Larvicides are not as effective as adulticides and, like adulticides, can prompt rapid evolution of resistance. However, Koella et al. (2009) suggest that this could actually benefit control, for the following reason. Resistance often comes with an evolutionary cost and, in Culex, resistant individuals may have shorter lifespans (Gazave et al. 2001). This is important because, as noted in §1, the adult mosquito lifespan is one of the most crucial factors affecting the transmission of malaria. Other evolutionary costs of resistance can include longer developmental times and smaller size as adults (which in turn is known to be correlated to longevity and biting rate; Koella et al. 2002; Lehmann et al. 2006). A combination of larviciding and the use of an LLA insecticide could lead to some particularly interesting future work owing to resistance to the two insecticides involving opposite changes in the traits associated with resistance (Koella et al. 2009). Costs of resistance to an LLA insecticide could, as emphasized by Read et al. (2009), lead to such an insecticide becoming completely evolution proof if resistance costs outweigh resistance benefits. The reason this becomes a realistic possibility is that the fitness gains of resistance to an LLA insecticide benefit only a few mosquitoes (those that live to old age), whereas the fitness costs (which may be additional mortality or reduced fecundity) are paid by all.

Finally, it should be emphasized that the strategy of killing only old, potentially infectious adult mosquitoes emphasizes disease control, and not insect control. As currently used, chemical insecticides control mosquitoes by killing individuals of all ages. Since mosquitoes are generally perceived as a nuisance, this is popular with the public. However, killing mosquitoes of all ages does increase the selection pressure for insecticide resistance (Read & Thomas 2009). There is thus a tradeoff between effective prevention of malaria transmission by mosquitoes, and having to live with mosquito bites involving no malarial transmission.


We would like to thank the Fields Institute for their kind support and hospitality during the Summer 2010 Thematic Programme on the Mathematics of Drug Resistance in Infectious Diseases. We also thank the referees for valuable comments.

  • Received August 3, 2010.
  • Accepted December 16, 2010.


View Abstract