## Abstract

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 §2*a*, the local and global stability of the equilibrium in which the resistant strain is established and the vulnerable strain is extinct. Then, in §2*b*, 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 §2*b*, 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
2.1where the density *v*(*t*,*a*) is defined such that for infinitesimal d*a* the number of vulnerables aged between *a* and *a*+d*a* is *v*(*t*,*a*) d*a*. 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,
2.2The total number of adult vulnerables is
and, using equation (2.2) and assuming that , it satisfies
2.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
2.4Thus, , and so equation (2.3) becomes the *V* equation of the system
2.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 . 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 §2*b*, 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
2.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τi}*b′(R***)>0. Then, the equilibrium (V,R)=(0,R***) of system (2.5) is locally asymptotically stable.*

### Proof.

Linearization of system (2.5) at the equilibrium (0,*R**) yields the partially decoupled system
2.7where . In view of equation (2.6), the first of these equations can be put in the form
2.8and, since *δ*>0, Kuang (1993, theorem 3.2.1) yields that as . The equation then becomes
Since *μ*_{a}>e^{−μiτi}*b*′(*R**)>0, as (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:
2.9The Nicholson’s blowflies birthrate, *b*(*M*)=*pM*e^{−qM}, 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].*

### Proof.

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),
2.10so that . A graphical argument using the properties of *b*(⋅) shows that e^{−μiτi}*b*_{max}/*μ*_{a}<*R*_{max}. Therefore, *M*(*t*)<*R*_{max} 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
We 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],**R**^{2}), but it remains valid when applied to the region 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τi}*b*(*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
where ≤ 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 (*V*_{t}(*θ*),*R*_{t}(*θ*)), *θ*∈[−*τ*_{i},0], is strictly order preserving with respect to <_{K}. Here, < means ≤ and ≠, and *V*_{t}(⋅) is the function *V*_{t}(*θ*)=*V* (*t*+*θ*), *θ*∈[−*τ*_{i},0]. So we need to check that implies that , i.e. that
with or on [−*τ*_{i},0], implies
with or , *θ*∈[−*τ*_{i},0]. This amounts to showing strict order preserving for the variable (*V*,−*R*), in the sense of the usual partial ordering in **R**^{2}. Letting and , we find the evolution equations for to be
2.11where
Routine differentiation, using that *b* is increasing and that , shows that *F*_{1} and *F*_{2} 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, and for *θ*∈[−*τ*_{i},0]. Then, the two solutions (*V*,*R*) and of system (2.5) agree on an interval of length *τ*_{i} and, by uniqueness, for all *t*>*t**. We show that this contradicts or on [−*τ*_{i},0] by showing that it implies and on [−*τ*_{i},0]. Introducing and , these functions are zero for *t*≥*Nτ*_{i}, where *N* is an integer such that *Nτ*_{i}≥*t**. Hence, from system (2.5),
for *t*≥(*N*−1)*τ*_{i}, where
Under the hypotheses, *G*_{1} is increasing with respect to *V* and decreasing with respect to *R*, while *G*_{2} is increasing with respect to *R* and decreasing with respect to *V* . Therefore, for *t*≥(*N*−1)*τ*_{i},
for some function . Since ∂*G*_{1}/∂*V* >0, this implies *ξ*(*t*)≤0, and so *ξ*(*t*)=0. We have shown that *ξ*(*t*)=0 for *t*≥*Nτ*_{i} 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 *G*_{2}. So, and 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τi}*b*′(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*)=−*μ*_{a}*R*(*t*)+e^{−μiτi}*b*(*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 and either (*V* (⋅),*R*(⋅)) or belongs to int (*X*^{+}×*X*^{+}), then 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
yields

The requirement that (*V* (⋅),*R*(⋅)) or belongs to int (*X*^{+}×*X*^{+}) guarantees that *φ*(*t*−*τ*_{i}) and cannot both be zero, so that . 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*)=−*μ*_{a}*R*(*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
which results from the use of the trial solution 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, *V*_{a} is the number of adult vulnerable mosquitoes and *V*_{o} is the number of old vulnerable mosquitoes. We propose the following model, which we call the LLA insecticide model:
2.12where *M*_{a}(*t*)=*V*_{a}(*t*)+*R*_{a}(*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 *V*_{a} 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
2.13For vulnerable mosquitoes at all stages of life, their age density *v*(*t*,*a*) satisfies
2.14where
2.15Differentiating equation (2.13) and using equations (2.14) and (2.15) gives
2.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

With *μ*(*a*) given by expression (2.15), this gives
2.17and
2.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*(*M*_{a}(*t*)+*M*_{o}(*t*)). As in the previous section, we deduce that
2.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 *V*_{a} equation of system (2.12).

The *V*_{a} equation of system (2.12) essentially formulates the change rate of the *V*_{a} 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
if this matching condition is met initially. So system (2.12) can be replaced by the differential equations for *V*_{o} and *R*_{o}, coupled with the above integral equation and its analogue for *R*_{a}. 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 §2*a*, 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, *V*_{a}=*V*_{o}=0 and , *R*_{o}=*R**_{o} where
2.20Thus, and we have a single equation for *R**_{o},
2.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
2.22For the linear stability of the equilibrium of system (2.12), we shall need the following additional condition:
2.23We now investigate this linear stability. Near the equilibrium , the linearization of the *V* equations is
2.24Setting and yields the linearization of the *R* equations,
2.25and
2.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
2.27We solve equation (2.27) together with the second of system (2.24) using the usual ansatz , where *c*_{1} and *c*_{2} are constants. After some algebra, this leads to the following useful form for the characteristic equation to be solved for *λ*:
2.28where
2.29In theorem 2.3, we prove that all roots of equation (2.28) have a negative real part, and this means that 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
2.30the first of which can be rewritten as
2.31This leads to the second characteristic equation
2.32where
2.33

We prove the following result on the linear stability of .

### Theorem 2.3

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

### Proof.

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
2.34Writing , it is necessary to show that, for *x*≥0,
2.35or, equivalently, that
2.36where . Let us first prove that, for a fixed ,
2.37Since , 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
so that
To 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 for *ξ*>0, so that *j*(*ξ*)>0 when *ξ*>0. Having established equation (2.37), to show (2.36) it is necessary to show that
2.38and this follows from the fact that the function is decreasing (see Gourley *et al.* 2007). We have shown inequality (2.34).

Now suppose there exists a root of equation (2.28) such that . Then,
2.39Writing *f*(*λ*) in the form
2.40where *g*(*x*)=(1−e^{−x})/*x*, we have, since ,
It is our claim that, for ,
2.41The right-hand side of this simplifies to because of (2.20). Therefore, in view of the inequality |*z*_{1}−*z*_{2}|≥|*z*_{1}|−|*z*_{2}| for *z*_{1},*z*_{2}∈**C**, it suffices to show that
2.42But, using the first equation of (2.20) and the definition of *g*, we can show that inequality (2.34) is equivalent to
and inequality (2.42) follows immediately, since . Using inequality (2.41) and the simplification of its right-hand side, we may now deduce from equation (2.39) that
using the second equation of (2.20). But for , this is impossible because it implies that 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
2.43and *f*_{1}(*λ*) is decreasing for *λ*≥0. Note that inequality (2.23) implies
so that
2.44Since inequalities (2.23) and (2.44) hold, we have
so inequality (2.43) holds. To show that *f*_{1}(*λ*) 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, *f*_{1}(*λ*) has a decreasing numerator and increasing denominator (note that , from inequality (2.23)). It follows that *f*_{1}(*λ*) is decreasing if
when *λ*≥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 μ*_{o}*=μ*_{a} *and δ>0. Then, the convergence of solutions to the equilibrium* *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
*2.45

### Proof.

The eigenvalue *λ*_{all} satisfies
2.46and, since *μ*_{o}=*μ*_{a}, *λ*_{old} satisfies equation (2.28), i.e.
2.47Define *f*_{1}(*λ*)=*μ*_{a}e^{−λτi} and
with *f*(*λ*) given by (2.29). Sketches of the curves , and and continuity arguments show that, if
2.48then inequality (2.45) follows. Since we assume *μ*_{o}=*μ*_{a}, we may add the equations in (2.20) to get , and so
Also, *f*_{1}(*λ*_{all})=*μ*_{a}e^{−λallτi}. So, to show that *f*_{2}(*λ*_{all})>*f*_{1}(*λ*_{all}), it is necessary to show that
which is equivalent to
and the truth of this follows from equation (2.46). It is easily checked that *f*_{2}(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*)=*pM*e^{−qM}, 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 1*a* 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. 1*b*), 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 *R*_{a}+*R*_{o} 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 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}.

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.

## Acknowledgements

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.

- This journal is © 2011 The Royal Society