A model of insect control with imperfect treatment

Insecticide spraying of housing units is an important control measure for vector-borne infections such as Chagas disease. However, some vectors may survive treatment, due to imperfect spraying by the operator or because they hide deep in the cracks or other places, and re-emerge in the same unit when the effect of the insecticide wears off. While several mathematical models of this phenomenon have been previously described and studied in the literature, the model presented here is more basic than existing ones. Thus it is more amenable to mathematical analysis, which is carried out here. In particular, we demonstrate that an initially very high spraying rate may push the system into a region of the state space with low endemic levels of infestation that can be maintained in the long run at relatively moderate cost, while in the absence of an aggressive initial intervention the same average cost would only allow a much less significant reduction in long-term infestation levels.


Introduction
Models of insecticide spraying often assume that treatment of a given unit will result in temporary removal of all insects from this unit. For example, this assumption was made in our previous paper [15], but such a scenario will not always occur. Here we relax this assumption and study a more realistic model that assumes that in some cases insects may survive treatment. After hiding deep in the cracks or similar places, these insects will reemerge after the effect of the insecticide has worn off.
In the literature on Chagas disease, re-emergence of triatomine infestation after treatment has been well documented. This phenomenon can be caused by imperfect spraying or because the insects could hide in bedding, cardboard boxes storing clothing, agricultural produce stored within the house, or other spaces that could not be sprayed or not appropriately exposed to sunlight when taken outside the house during insecticide application [8].
The model presented here could apply to any of the causes for re-emergence described in the previous paragraph, except insecticide resistance (see our discussion in Section 4). While inspired by our interest in controlling the spread of Chagas disease, the model and our results would apply to a wide variety of systems where housing units are (re)-infested by insects. Insecticide treatment with less than perfect effectiveness has already been studied in several mathematical and simulation models [1,3,16,[23][24][25]. However, our imperfect treatment model is more basic than previously published ones in the sense that it has very few variables and parameters and is entirely based on ODEs, as opposed, to, for example, delay differential equations. Thus it is more amenable to analytic exploration.
One of the main goals of our modelling is assessment of cost-effectiveness of insecticide treatment strategies that eventually could become a basis for public policy recommendations. Our major finding is that the dual-rate effect that was shown to occur in the models studied in [15] also does occur in the imperfect treatment model presented here. This has important implications for control efforts, as in the presence of this effect a strategy of adopting a highly aggressive short-term intervention at the beginning of a sustained insecticide control programme may be more cost-effective in the long run than spraying at a fixed rate without a heightened initial effort.
The paper is organized as follows: In Section 2 we define the imperfect treatment model that is the focus of the present paper and study its basic mathematical properties. In Section 3 we prove that the above-mentioned dual-rate effect can occur in the imperfect treatment model as long as the parameters satisfy certain conditions. Section 4 contains a summary of our findings, discussion of their significance, and open problems for follow-up research.

The imperfect treatment model
Hosts in this model are units that can be infested by insects. A unit could be a house, a peridomestic area, or the house together with the peridomestic area. A similar approach was used in [2,15,26], and to a certain degree in [12].
We do not make a distinction between insects that carry and that do not carry pathogens. In a sense, we are simply treating the presence of insects as posing some risk of disease transmission and are interested in strategies for diminishing this risk by controlling the overall insect population. A unit is considered infested if it contains a viable colony of insects, which should be understood as a colony that will not die out unless the unit is treated with insecticide. In field studies, detection of at least one live insect within the unit is usually considered a good proxy for presence of a viable colony (see, for example, [7,15] for more discussion of this issue).
An uninfested unit can become infested by migration of insects from an infested unit or from a reservoir such as a sylvatic area, unless the unit has been recently treated with insecticide and is temporarily immune to (re)infestation. The effectiveness of the insecticide is assumed to decay exponentially. Only infested units are assumed to receive insecticide treatment. Moreover, the imperfect treatment model allows for the possibility that some insects may survive treatment and temporarily hide within the treated unit while  The rate at which effectiveness of the insecticide wears off r Spraying rate of units α The fraction of infested units that get successfully treated the insecticide remains effective. We conceptualize the model as an SIRS model with four compartments: uninfested (S), infested (I), and two compartments R 1 and R 2 of temporarily removed units. Here R 2 contains successfully treated and temporarily immunized units that will enter the S-compartment when the effect of the insecticide wears off. In contrast, compartment R 1 contains units for which treatments was only partially successful, so that viable colonies survive inside these units. While in R 1 , units are assumed to be protected from infestation from other sources and not themselves to be sources of infestation of other units. However, once the effect of the insecticide falls below a certain threshold, units in R 1 will move into the infested compartment I. Figure 1 illustrates this description, and also shows the rates of movement between these compartments. Table 1 lists all parameters and their meanings. The imperfect treatment model itself is an ODE-based mass-action model without demographics (building of new units and destruction/abandoning of old ones). The variables represent numbers (rather than proportions) of units in the compartments. The total number of units will be denoted by m. As it seems plausible that the effect of the insecticide wears off at the same rate w in the R 1 and R 2 compartments, we defined and analysed the model only under this assumption. It may be of some mathematical interest to explore in the future a generalization of the model where these rates can be different for compartments R 1 and R 2 .
Let us remark that while our model is not about disease transmission and the units in our model are not organisms, we will freely use the terminology of mathematical epidemiology, such as 'basic reproductive ratio'. Similar to ODE models in mathematical epidemiology, our model can only give an approximation to the stochastic dynamics of numbers that in reality are integers. This approximation will be reasonable accurate if the population size m represents units in large batches. In our case, the underlying stochastic process would involve detection of insects in a unit and treatment thereafter at random times after (re)infestation or re-emergence. Note that this is different from the assumption of insecticide treatment at fixed time intervals that has been studied, for example, in [1,23,24] and also in one of the models of [14].
The differential equations describing the dynamics of the population in each compartment are: The basic model that was investigated in [15] can be considered a special case of model (1) with α = 1 for the restriction to states with R 1 = 0.
Note that S = m − I − R 1 − R 2 . Thus we can reduce (1) to the following threedimensional model: The biologically feasible region of model (2) is given by To show that the region is forward invariant, consider the derivatives on its boundary: Direct inspection of the system (2) shows that dI/dt| i=0 ≥ 0. Similarly, dR 1 Thus when P(0) = I(0) + R 1 (0) + R 2 (0) = m, then dP/dt = −wR 2 ≤ 0. For the analysis of the dynamics of model (2), first notice that along all trajectories the ratio between R 1 and R 2 approaches a fixed value when t → ∞. Specifically, let

Corollary 2.2: (a) The models (2) do not have nonconstant periodic orbits.
(b) Let − ⊆ be some closed forward invariant region such that − ∩ = is simply connected. If for some choice of parameters the model (2) has a unique equilibrium in − that is locally asymptotically stable, then this equilibrium must be globally asymptotically stable on − .
Proof: For the proof of part (a), see Appendix 4.
For the proof of part (b), notice that by Proposition 2.1, the full system (2) and its restriction to = must have the same equilibria. Since the set = is simply connected, the analogue of part (b) for the restriction of the dynamics to = follow from the Poincaré-Bendixson Theorem. Now part (b) is implied by Proposition 2.1 and geometric constraints on the trajectories.
Since for α = 1 model (2) reduces to a reparametrization of the basic model of [15] for which the results of this paper were already proved, we can focus here on the case 0 ≤ α < 1. Then on = we have R 2 = (α/(1 − α))R 1 and (2) reduces to a two-dimensional system The following theorem summarizes the properties of our model when there is no invasion from the sylvatic areas, that is, when c = 0. (2) and assume c = 0. Then (i) When βm − αr < 0, then the infestation-free equilibrium IFE = (I * , R * 1 , R * 2 ) = (0, 0, 0) is the only equilibrium in the biologically feasible region . It is both locally and globally asymptotically stable.

Theorem 2.3: Consider the model
(ii) When 0 < βm − αr, then the IFE becomes unstable. In addition to the IFE there exists a unique endemic equilibrium which is locally asymptotically stable. The EE will be asymptotically approached by all trajectories that start with I(0) + R 1 (0) > 0.

Proof:
The expression for the IFE follows from our assumption that w > 0 so that when I * = 0 all units will eventually revert to the susceptible state and we must have By Proposition 2.1, the equilibria of the models (2) and (3) coincide and have the same stability properties. In order to find endemic equilibria of both models, we set the second line of (3) to zero, which gives After substituting this into the first line of (3) and simplifying the resulting expression we obtain This shows that an endemic equilibrium must satisfy I * = w(βm − αr)/β(w + r). In particular, and endemic equilibrium exists in the biologically feasible region if, and only if, 0 < βm − αr. It will be unique in this case, and its coordinates must satisfy (4) by the definition of = .
Local stability of the equilibria can be determined by analysing the characteristic polynomials of the Jacobians. The Jacobian matrix of (3) is For c = 0 the Jacobian matrix at the IFE where I = R 1 = 0 reduces to The eigenvalues are At least one of the roots of (6) has a positive real part if, and only if, 0 < βm − αr. Hence the infestation-free equilibrium is locally asymptotically stable when βm − αr < 0. It is unstable whenever 0 < βm − αr.
Similarly, for c = 0, the Jacobian matrix at EE where R * After substituting the value of I * from (4) and algebraic simplifications we find that the characteristic equation of the Jacobian matrix is All the coefficients of (7) are positive when 0 < βm − αr. It follows from the quadratic equation that the endemic equilibrium is locally asymptotically stable whenever it is biologically feasible. Global stability of the IFE for spraying rates r > βm/α follows from part (i) and Corollary 2.2(b).
To prove global stability of the EE for initial conditions with I(0) > 0 and small spraying rates with αr < βm, we cannot rely on local stability and Corollary 2.2(b) only, since the IFE is not repelling. In our argument we again restrict ourselves to the case when 1 − α > 0, as the case α = 1 is already covered by the analogous result for the basic model in [15].
For ε ≥ 0 define the set IR≥ε and its boundary IR=ε as Choose ε, δ > 0 such that By Proposition 2.1, the set Since we assume c = 0 here, on the boundary IR=ε ∩ δ we get from (2) and (8): We conclude that the set IR≥ε ∩ δ is a closed forward invariant subset of that does not contain the IFE, but will contain the EE if ε is sufficiently small. The set IR≥ε ∩ δ ∩ = is also simply connected. Since the set IR=0 is invariant, for every initial condition outside IR=0 we can find ε, δ > 0 that satisfy (8) such that its trajectory eventually enters IR≥ε ∩ δ . Now the result follows from Corollary 2.2(b). This makes intuitive sense, as α = 0 implies that each spraying event leaves a viable colony in the treated unit that will eventually re-emerge.
(b) The basic reproductive ratio in this model is R 0 = βm/αr. A calculation based on the next-generation matrix is given Appendix 4. Note that additional secondary infestations may be caused by an already treated unit if that treatment was unsuccessful. As a consequence, when βm/r < 1 < R 0 , the infestation level in a population that consists entirely of susceptible and a small fraction of infested units will initially decrease before approaching the value I * of the EE from below. This makes it difficult to assess whether initial eradication efforts are sufficient to drive the population towards the IFE. (c) Taken together, parts (i) and (ii) of Theorem 2.3 imply that the system undergoes a bifurcation when r = βm/α, or, equivalently, at R 0 = 1. We have analytically confirmed that this is a transcritical forward bifurcation.
The following result summarizes the properties of the model when there is invasion from the sylvatic areas, that is, when c > 0. (2) and assume c > 0. Then (i) There exists a unique equilibrium EE = (I * , R * 1 , R * 2 ) in . It is given by

Theorem 2.4: Consider the model
(ii) The equilibrium EE is endemic and both locally and globally asymptotically stable on .
Proof: Again, we work with the restriction (3) of the model to = and rely on Proposition 2.1 for transfer of our results to the model (2). At equilibrium of (3), we have By the second line of (10) After substituting this into the first line of (10) and then performing straightforward algebraic transformations we obtain Since β(w + r) = 0, it follows that The only biologically feasible root of the resulting quadratic polynomial is The expressions for R * 1 and R * 2 then follow by setting the the second and the third line of (2) to 0. To show that I * ∈ (0, m), let Then f (0) = wcm > 0 and f (m) = −βrm 2 − αrwm − rcm < 0. Thus I * ∈ (0, m), and EE must be an endemic equilibrium in the interior of = , and it must be unique. This proves part (i).
Notice that the expression (4) for the EE when it exists in the model without invasion from sylvatic areas can be obtained by substituting c = 0 into (9).
For part (ii), it suffices to prove local asymptotic stability of the EE. Since there are no other equilibria when c > 0, global asymptotic stability then follows from Proposition 2.1 and Corollary 2.2(b).
For c > 0, the Jacobian matrix of (3) at EE becomes Its characteristic equation can be expressed as Details are given in Appendix 4. Notice that when βm − αr > 0, then (7) becomes the special case of these expressions for c = 0. The characteristic polynomial is quadratic and has positive coefficients. So all the roots have negative real parts, and local asymptotic stability of the EE follows.

Remark 2.2:
Note that from the first line of Equation (9) we get Since the EE is unique, this implies that for any given positive threshold there exists a sufficiently high spraying rate that will drive the infestation level at the EE below this threshold.

A dual-rate effect
Consider the model introduced in the previous section and a spraying rate r. The results of Section 2 indicate that when r is kept constant, in the long run the numbers of susceptible, infested, and removed units will approach an equilibrium, with I * (r) denoting the equilibrium number of infested units. Thus we can conceptualize the long-term cost of insecticide treatment in these models as follows: In particular, when c = 0 and if the spraying rate is sufficiently high to drive the system to the IFE with I * = 0, the long-term cost will be considered 0. It is intuitively clear that the equilibrium value I * (r) is non-increasing with respect to r. We have also confirmed this analytically (details not included here).
However, under certain conditions the dependence of C on r is not monotone.

Definition 3.1 (Dual-rate effect):
Let I * (r) be the numbers of infested units at endemic equilibrium and let C(r) be defined as above. A dual-rate effect occurs if there exist two different spraying rates r 1 and r 2 such that C(r 1 ) = C(r 2 ) and I * (r 1 ) < I * (r 2 ).
We focus here on endemic equilibria, even when an IFE exists and can be maintained at no cost, as will be the case for c = 0 by Theorem 2.3. The dual-rate effect may then still be of relevance if the EE exists and there is some upper bound on a biologically feasible spraying rate that may prevent approach to the IFE.
We will use the following observation.

Proposition 3.2:
If the cost function C(r) has a strict local maximum at some r ∈ (0, ∞), then the dual-rate effect occurs.
An example of numerical simulations is displayed in Figure 2. The curves are only slightly different and qualitatively the same as the ones for the basic model shown in Figure 2 of [15] for similar parameter settings. They are consistent with the theoretical predictions of the next two theorems, which are also analogues to results about the basic model that were proved in [15].

Theorem 3.4:
Consider the model (2). For any fixed parameters α, β, c, w > 0, the dual-rate effect occurs whenever m is sufficiently large. (9) and (13), the cost at equilibrium is

Proof: By Equations
For r = 0 we have We want to find a second point r = r 1 such that C (r 1 ) < 0. Choose Then r 1 > 0 as long as βm > c, and when m is sufficiently large (see (A5) and (A6) of Appendix 4 for details), Hence C(r) has a strict local maximum in the interval [0, r 1 ]. The conclusion now follows from Proposition 3.2.

Remark 3.1:
In order for m to work for a given choice of the other parameters, our proof requires that m > c/β and another very technical condition that can be deduced from (A5) and (A6) at the end of Appendix 4. These conditions could presumably be weakened. However, some assumptions on the parameters are needed in Theorem 3.4. This already follows from the results in [15] on the basic model, which is a special case of the imperfect treatment model. It is also illustrated in Figure 2, where the effect is shown to occur when the ratio c/β is small relative to m, but not when the ratio c/β is large relative to m.
It is of interest to compare the strength of the dual-rate effect for the perfect treatment model of [15] (when α = 1) and the imperfect treatment model in this paper when α < 1 and all other parameters of the models are same. Figure 3 illustrates the difference for α = 0.4. It can be seen that with imperfect treatment the peak that represents the least optimal spraying rate is located further to the right. We have confirmed this analytically (details not included). Moreover, with imperfect treatment the curve is flatter than with perfect treatment. This means that with imperfect treatment, the effect may either not permit the same magnitude of cost-savings as in the basic model (see right panel where c > 0) or will require substantially larger spraying rates (see the left panel where c = 0).
In this context, our theorems show that while the dual-rate effect may be weaker with imperfect treatment, it nevertheless persists in some form for suitable choices of the other parameters whenever α > 0, that is, when there is a positive probability that treatment will eliminate the colony in a given unit. The assumption that α > 0 is essential here, as the following result shows:

Summary, discussion, and open problems for future research
Insecticide spraying remains an important control measure for vector-borne diseases. Here and in [15] we investigated ODE-based models of type SIRS with a reservoir. These models reflect important aspects of (re)infestation by insects if the effectiveness of the insecticide wears off over time, the insecticide treatment is applied to units where insects can establish colonies and/or transmit pathogens to humans, and infestation can originate both from other units and a reservoir, such as sylvatic areas. All three of these aspects play important roles in the ecology of triatomines, the insect vectors for the parasites that cause Chagas disease. In contrast to the models of [15], the model studied here also accounts for the possibility that some insects may survive treatment and re-emerge in the same unit after the effectiveness of the insecticide has worn off. Such scenarios have been previously documented in the the empirical literature [8]. Less than perfect effectiveness of the insecticide has already been studied in several mathematical and simulation models [1,3,16,[23][24][25]. However, our imperfect treatment model is more basic than previously published ones and thus more amenable to analytic exploration.
In particular, we were able to derive explicit formulas for the equilibria in our models and characterize their stability. When there is no invasion from the (sylvatic) reservoir (c = 0), then a sufficiently high spraying rate guarantees that the system will be driven to the infestation-free equilibrium IFE. When c > 0, then there is always a unique equilibrium EE that is endemic and globally asymptotically stable. For any given positive threshold there exists a sufficiently high spraying rate that will drive the infestation level I * at the EE below this threshold.
Consider a situation where a community has decided to initiate an insecticide control programme that would in the long run keep the equilibrium infestation level below some threshold I * < I . The threshold I should be as low as feasible within the available longterm budget. Would it then be possible to maintain a lower threshold in the long run, within the same budget, by using initially very high spraying rates that would not be sustainable in the long run? In essence, could such an initial short-term and very aggressive intervention drive the system into the region where maintaining a lower threshold is feasible in the long run than would be feasible without such an initial effort? This question had not been addressed in the literature prior to our work (see [15] for a comprehensive discussion). The answer will be positive whenever the dual-rate effect of Definition 3.1 occurs in the model. Our results show that this effect always occurs for any fixed choice of positive values for the other parameters when the number of units m is sufficiently large. This finding is analogous to the results of [15] for models with perfect treatment.
If the effect does occur, then one would like to recommend initially aggressive spraying to lower long-term cost. As was already pointed out in [15], some caution is needed when making such a recommendation, because the effect depends on the parameters and assumptions of our model; see Remark 3.1. There is an additional caveat here: One would need to taylor the initial intervention so that it drives the total number of units that carry insects below a certain threshold. But this would need to be a threshold for I + R 1 . Thus it will in general be different from the target level I , and it may not be directly observable. Moreover, the model parameters may be more difficult to estimate than when treatment is always perfect; see, in particular, Remark 2.1(b).
We constructed our imperfect treatment model to make it as simple as possible. This very simplicity is an advantage insofar as it made the model amenable to analytic study. On the other hand, it opens up a number of ways in which the model can be elaborated and made more biologically realistic. From our point of view the most interesting question would be whether the dual-rate effect would still occur in such more elaborate versions.
Some promising directions for such elaborations apply both the imperfect treatment model and the models studied in [15]. They include incorporating of other types of control measures in addition to insecticide spraying, of demographics (constructing new housing units and abandoning old ones), of various nonlinearities and/or a discount factor into the cost function (13), consideration of spraying at fixed intervals rather than continuously (as in [1,23,24] and one model of [14]), relaxing the assumption of exponential decay of the insecticide at a fixed rate, and explicitly incorporating a cost of insect detection. For a more detailed discussion of these options see [15].
Let us focus here on two aspects that apply more specifically to possible refinements of the imperfect treatment model.
One important and interesting question is how to adapt our model to the study of situations where unsuccessful treatment is caused, at least in part, by the development of insecticide resistance. This scenario has been empirically documented in the literature [9,10,13,[19][20][21][22][27][28][29]. In this case our model may not apply in its present form, for at least three reasons: First, (partial) insecticide resistance of the colony within a given unit would in effect lower the constant α for that unit. Thus unless such (partial) resistance is universal for all infested units or the only reason for unsuccessful treatment, we might need a model with more compartments and rate constants α 1 , α 2 . Second, a resistant colony may not undergo a period when insects hide and don't infest other units, or such a period might be shortened relative to units with non-resistant strains. The latter can be dealt with by the generalization of our model to the case where the rates at which units leave the compartments R 1 and R 2 are different, so that we have two rate constants w 1 > w 2 instead of one constant w. But this would not all by itself work for modelling (partial) insecticide resistance without adding compartments to our model. Third, due to strong selection pressure, insecticide resistance of a colony that survives repeated spraying events is likely to increase over time, so that the coefficients α 1 , w 1 would need to vary over time as well.
The previous paragraph outlines a model, but makes an implicit assumption that the same insecticide is always used. When the community switches to a new spraying agent, the situation goes back to one that conforms to the imperfect treatment model studied here, at least for a while. This brings us to our second open problem: How to incorporate behaviour changes, such as switching to a new spraying agent or taking protective measures, like sealing cracks or installing insect screens, that would lower a unit's susceptibility, into the model? Such behaviour changes would be triggered by awareness that may result either from direct experience of previous infestations or from spreading by word of mouth between owners of the units. Models with spread of awareness have been studied for infectious diseases of type SIS in [11]. It would be of great interest to extend this work to models where the underlying dynamics is of type SIRS with a reservoir.