In AndySouth/resistance: an insecticide resistance population genetics model with 2 loci and 2 insecticides

Andy South * and Ian M. Hastings

Department of Parasitology, Liverpool School of Tropical Medicine, Liverpool L3 5QA, UK.

southandy@gmail.com

ian.hastings@lstmed.ac.uk

* Corresponding author

\pagebreak

#!!!! PAPER PROTOCOL
#!!!! submission MS, knit to word_document 
#!!!!  in WORD: delete figs, set double space, add line numbers 
#!!!!  (spacing & line numbers work for pdf output but not for word)
library(knitr)
opts_chunk$set(dev="tiff",
               dev.args=list(compression="lzw"),
               dpi=300,
               cache=FALSE,
               fig.path='pap2figs/')

  library(resistance)
  library(ggplot2)

Abstract

Background

Insecticide resistance threatens effective vector control, especially for mosquitoes and malaria. To manage resistance, recommended insecticide use strategies include mixtures, sequences and rotations. New insecticides are being developed and there is an opportunity to develop use strategies that limit the evolution of further resistance in the short term. A 2013 review of modelling and empirical studies of resistance points to the advantages of mixtures. However, there is limited recent, accessible modelling work addressing the evolution of resistance under different operational strategies. We see an opportunity to improve the level of mechanistic understanding within the operational community of how insecticide resistance can be expected to evolve in response to different strategies.

We provide a concise, accessible description of a flexible model of the evolution of insecticide resistance. We use the model to develop a mechanistic picture of the evolution of insecticide resistance and how it is likely to respond to potential insecticide use strategies. The aim is to reach an audience unlikely to read a more detailed modelling paper. The model itself, as described here, represents two independent genes coding for resistance to two insecticides. This allows us to represent the use of insecticides in isolation, sequence and mixtures.

Results

We use the model to demonstrate the evolution of resistance under different scenarios and how this fits with intuitive reasoning about selection pressure. We show that using an insecticide in a mixture, relative to alone, always prompts slower evolution of resistance to that insecticide. However, when resistance to both insecticides is considered, resistance thresholds may be reached later for a sequence relative to a mixture. Increasing the ability of insecticides to kill susceptible mosquitoes (effectiveness), has the most influence on favouring a mixture over a sequence because one highly effective insecticide provides more protection to another in a mixture.

Conclusions

Our model provides an accessible description of the process of insecticide resistance evolution and how it is likely to respond to insecticide use. We provide a simple user-interface allowing further exploration (https://andysouth.shinyapps.io/resistmob2). These tools can contribute to an improved discussion about operational decisions in insecticide resistance management.

Keywords

insecticide resistance; public health; mosquitoes; vector-borne diseases; infectious diseases; malaria; population genetics

Background

Insecticide resistance is a problem for malaria control[@WHO2012][@Ranson2016][@Hemingway2016] other vector borne diseases [@IRAC2011] and agriculture [@FAO2012]. Malaria alone still results in hundreds of thousands of deaths per year. Recent malaria control efforts have centred on treated bed nets and indoor residual spraying, both reliant on insecticides. Treated nets were recently estimated to contribute 68% and indoor residual spraying 13% to averting more than 500 million falciparum malaria cases between 2000 and 2015 [@Bhatt2015]. A recent malaria transmission model [@Churcher2016] predicts that even low pyrethroid resistance is likely to increase malaria incidence in Africa by reducing the performance of bed nets.

The WHO produced a Global Plan for Insecticide Resistance Management in malaria vectors (GPIRM)[@WHO2012] which includes recommendations on operational strategies for managing resistance including the use of insecticide mixtures when they become available. Efforts are under way to develop new insecticides that will be effective in the light of existing resistance and allow additional options within insecticide resistance management. The Innovative Vector Control Consortium (IVCC) was set up in 2005 to work with industry to develop new vector control tools and particularly new insecticides with new modes of action to address insecticide resistance in disease transmitting mosquitoes [@Hemingway2006][@IVCC2016]. Three new insecticides are now in development [@IVCC2016] and likely to be available after 2020 [@Ranson2016]. It is important that decisions about how best to use the new insecticides, with a clear mindset of delaying the evolution of resistance, are made before insecticides are released [@Hemingway2016].

A recent comprehensive review of strategies to avoid resistance evolution across pesticides and drugs [@Consortium2013] concluded that mixtures (combination of molecules) are usually the best resistance management strategy. This was based on both empirical and modelling work. Modelling studies have investigated the evolution of insecticide resistance in insecticide mixtures including in a public health context e.g. [@Curtis1985][@Mani1985][@Roush1989] but much of the work was done more than 25 years ago and there remained some confusion about the results [@Levick2017]. In an earlier paper [@Levick2017] we described the technical details of a flexible model used to investigate the relative benefits of mixtures and sequences. For that paper we investigated thousands of scenarios to explore potential outcomes. In contrast here we provide an accessible summary of the model and use selected parameter values to describe mechanistically how the evolution of resistance can be expected to respond. This mechanistic understanding can contribute to the debate on the relative merits of different insecticide strategies and extend existing frameworks [@Roush1989][@IRAC2011][@FAO2012][@WHO2012][@Consortium2013]. By making the model more accessible we also hope to help bridge gaps between academia and policy in insecticide resistance management [@Huijben2017].

Our modelling approach focuses on the change in frequency of a single resistance gene for each insecticide. We assume that genes conferring resistance are already present in the population and that resistance to each insecticide is coded by a single gene. These assumptions are consistent with previous modelling work [@Consortium2013]. In support of the first assumption insecticide resistance has been termed 'pre-adaptive' meaning that resistance alleles are present in very low numbers, most likely as a mutation-selection balance, even prior to exposure to novel insecticides [@Ffrench-Constant2013][@Liu2015] as has been documented in blowflies [@Hartley2006]. In support of using a single gene representation of resistance evidence suggests that whilst polygenic resistance is common, control failures, particularly those due to target site resistance, are mostly due to single major genes [@Ffrench-Constant2013].

Methods

The simulation represents a population of randomly mixing individuals using standard population genetic approaches to avoid the need to follow every individual. One gene or locus is represented for each insecticide. Each individual has two copies of the gene and there are two potential alleles per locus. Each allele confers either resistance or susceptibility to the insecticide. Thus individuals can either have both resistance genes (RR homozygous), both susceptible genes (SS) or one of each (SR, heterozygous). The combination of the two genes is termed the genotype. There are 3 genotypes when considering one insecticide and 9 genotypes (e.g. SSRR) when considering two insecticides (Fig 1). Fitness is the main currency of the model, representing how much each genotype survives and reproduces. The fitness of the SS genotype in the absence of the insecticide is used as the reference genotype and assigned a value of 1. All other fitnesses are assumed to be less than 1, determined by alleles, exposure to the insecticide and other inputs that can be set in the model. These inputs include the effectiveness of the insecticide against SS genotypes, the dominance of the resistance allele and the ability of resistance to restore fitness of the RR genotype in the presence of the insecticide (see Table 1).

Table 1. Parameters influencing the development of insecticide resistance

Parameter | Description ---------------------- | ---------------------------------------------------- 1. Effectiveness | proportion of susceptible (SS) insects killed by exposure to insecticide 2. Exposure | proportion of insects exposed to insecticide 3. Resistance restoration | ability of resistance (RR) to restore fitness of insects exposed to insecticide 4. Dominance of resistance | sets fitness of heterozygous (SR) insects between that of SS & RR in presence of insecticide 5. Frequency | frequency of resistance alleles within the population 6. Cost of resistance | fitness of resistant (RR) insects in absence of insecticide 7. Dominance of cost | sets fitness of heterozygous (SR) insects between that of SS & RR in absence of insecticide

The calculation of fitness for each genotype when considering just a single insecticide is illustrated in Figure 2. In the model itself this is repeated for both insecticides and the fitnesses multiplied. Firstly the exposure input determines the proportion of the population in the left and right panels (exposed and not exposed). For those that are exposed (left panel) insecticide effectiveness sets the fitness for SS and resistance restoration 'restores' a portion of the fitness for RR. A value of 1 for resistance restoration would bring the fitness of RR back to 1, equivalent to the unexposed SS and meaning the RR genotype is completely unaffected by the insecticide. A value of 0 for resistance restoration would lead to the RR having the same fitness as SS in the presence of the insecticide and effectively there would be no resistance. We use 'Resistance restoration' to represent how much the RR genotype restores fitness in the presence of the insecticide. Cost of resistance determines the fitness of RR in the portion of the population not exposed to the insecticide.

Dominance of resistance determines how the fitness for SR sits between that of SS and RR for those that are exposed to the insecticide. For those that are not exposed, fitness of SS is set to 1 by definition and resistance cost determines the fitness of RR. Dominance of cost determines how the fitness for SR sits between that of SS and RR in the absence of the insecticide (Fig 2).

The simulation proceeds through non-overlapping generations. In each generation selection is represented by multiplying genotype frequencies in the population by their relative fitness, the latter depending on their absolute fitness (Fig 2) and level of exposure. This acts to make the fitter alleles more common over time. Separate sexes and standard sexual reproduction with recombination is included. There is also the potential to set exposure rates to be different for males and females to represent their different behaviours. Females are more likely to come into contact with insecticides on walls and nets when seeking to feed on humans compared to males that do not. We only explore that input briefly here because a previous analysis indicated it was not a major factor determining differences between insecticide-use strategies [@Levick2017].

So far this description just considers a single insecticide and associated resistance allele. In the model itself a second insecticide and resistance allele are included and overall fitness is calculated by multiplying the two results. This allows the model to represent populations exposed to two insecticides together in a mixture.

The model outputs the change in resistance allele frequencies over time measured in generations. We extract the number of generations at which a resistance allele frequency of 50% is reached and term this 'time-to-resistance'. Using a 50% threshold is consistent with previous modelling analyses e.g. [@Curtis1985][@Birget2015]. Preliminary investigations showed that using other potential thresholds of 25 or 75% gave very similar results.

We investigate three insecticide use strategies.

single insecticide
two insecticides used in sequence, replacing the first with a second once the frequency of resistance alleles reaches a threshold of 50%
two insecticides used in a mixture (concurrently)

A major aim of this paper is to use modern techniques to make the modelling accessible to a wider community. The model is implemented in R [@RCoreTeam2017], the code is open-source and hosted on Github [@South2017] including the code to generate both the figures and text of this paper relying on the packages knitr[@Xie2017], shiny[@Chang2017] and ggplot2[@Wickham2009]. We provide an on-line user interface [@South2017a] (Fig 3) accessing the same code and enabling readers with no coding experience to change inputs and run the model themselves.

Results

We start by describing the results of model runs with a single insecticide which are more straightforward and move on to describing results of using two insecticides in sequences and mixtures.

Single insecticide

For single insecticide use, higher values of insecticide effectiveness, exposure, resistance restoration and dominance of resistance all resulted in faster resistance spread (Fig 4 A-D) and thus fewer generations to reach a resistance threshold of 0.5. This makes intuitive sense as all increase the impact of the insecticide and thus increase selection for resistance. Similarly, higher values of starting resistance frequency also resulted in shorter times to resistance thresholds (Fig 5A) through the even simpler mechanism that a smaller change in frequency is required to reach the threshold.

Cost of resistance (Fig 5B) and its dominance (not shown) had the opposite effect. Higher values led to a slower development of resistance. Again this is intuitive because costs reduce the advantage of resistance and decrease selection. These effects of inputs on single insecticide use are summarised in Table 2.

Table 2. Effect of inputs on resistance when insecticides used singly or in sequence

Parameter to increase | effect on resistance evolution | Mechanism ----------------------- | ------ | --------------------------------------------- 1. Effectiveness | faster | reduced fitness of SS & SR exposed to insecticide 2. Exposure | faster | reduced fitness of SS & SR overall 3. Dominance of resistance | faster | increased fitness of SR exposed to insecticide 4. Resistance restoration | faster | increased fitness of RR & SR exposed to insecticide 5. Frequency | faster | less change needed to reach resistance threshold 6. Cost of resistance | slower | reduced fitness of RR & SR not exposed to insecticide

Two insecticides

When two insecticides are used we can compare the relative performance of concurrent use in a mixture with sequential use. Each panel in Figures 6-10 compares mixture and sequence for a single combination of inputs. To help understand the mechanisms influencing this relative performance we start with a single base scenario and demonstrate changing inputs individually. In this base scenario we set effectiveness, exposure, resistance restoration, and dominance at 0.5, the starting frequencies of resistance to 0.01 and cost of resistance to 0. Resistance arises slower for sequential than mixture for this base scenario (Fig 6A). Resistance to both insecticides in the mixture follow the same path and reach the threshold at the same time (as would be expected given that they have identical input parameters). When used in sequence the curve (dashed) for each insecticide individually is steeper than in the mixture but, because they happen one after the other, it takes longer for both to reach the resistance threshold (Fig 6A). Thus the time-to-resistance for the mixture divided by the sequence is less than 1 at 0.8.

Insecticide exposure and effectiveness

Increasing exposure (the proportion of insects that come into contact with the insecticides) (Fig 5C) decreases the time-to-resistance for both the sequential and mixture strategy. The effect on the mixture is greater so time-to-resistance remains longer for the sequence (the ratio decreasing from 0.8 to 0.6). In contrast keeping the exposure constant and increasing effectiveness (the proportion of SS insects that are killed by contact) of one of the insecticides (Fig 6B) results in a longer time-to-resistance for the mixture relative to the sequence (and the ratio increasing from 0.8 to 1.2). The longer time-to-resistance for the mixture results from a slower increase in resistance for the less effective insecticide. Resistance to the less effective insecticide in the mixture increases slowly initially and then speeds up after the more effective insecticide has reached the resistance threshold. This points to a mechanism whereby the more effective insecticide initially slows the rate of evolution of resistance to the less effective one.

The result is that a more effective insecticide increases time-to-resistance when used in a mixture (compare the solid lines in Fig 6A & B). This is opposite to what happens when used in sequence (compare the red dashed line in Fig 6A & B) and in isolation (Fig 4A), when more effective insecticides shortened times-to-resistance.

In summary increasing effectiveness can favour mixtures moving from column 1 to 2 in Figure 6, and increasing exposure can favour sequences moving from row 1 to 2 in Figure 6. We do not show the results here, but reduced male exposure with respect to female, as would be expected for male mosquitoes not seeking blood meals, had a similar effect to reducing overall exposure.

Increasing the effectiveness further of either or both insecticides increases the performance of mixtures over sequences in terms of time-to-resistance. Figure 7A uses the same inputs as Fig 6B. From this scenario a higher effectiveness for insecticide2 (Fig 7C) or insecticide 1 (Fig 7B) or both (Fig 7D) all result in a greater positive difference in time-to-resistance for mixture over sequence. This again points to how more effective insecticides can slow resistance to another when used in a mixture.

Resistance restoration and dominance

Increasing resistance restoration or its dominance decreases time-to-resistance for both sequences and mixtures (Fig 8) in the same way that it did for sole use (Fig 4). The result is that increasing either or both of dominance and resistance restoration (Figs 8B,C,D) does not change the relative time-to-resistance for mixtures and sequences from that in the base scenario (Fig 8A) where it takes resistance longer to develop for the sequence. So, unlike effectiveness and exposure, the levels of resistance restoration and dominance do not effect whether a mixture or a sequence are likely to be favoured.

Starting frequencies of resistance

Changing the starting frequency of resistance had a similar effect on time to resistance for both sequences and mixtures and thus had little effect on their relative performance. For example taking the base scenario and reducing the starting frequency of one resistance allele did not change from sequence being favoured (compare Fig 9C to 9A). Similarly taking a scenario in which time-to-resistance is longer for a mixture (Fig 9B) and decreasing the starting frequency of that resistance allele did not change the fact that the mixture was favoured (Fig 9D). In these analyses we set starting frequencies of resistance relatively high at 0.01 or 0.001. In our more detailed sensitivity analysis [@Levick2017] the lower limit for starting frequencies was 0.0001 and still it was not a key parameter in determining the relative performance of mixtures and sequences.

Costs of resistance

Increasing the cost of resistance to one insecticide in a mixture leads resistance to that insecticide to evolve more slowly (compare the solid red line to the solid blue line in Figure 10C, the red line is for the resistance with the higher cost). The resistance with the higher cost also increases more slowly in the mixture relative to when it is used in sequence (compare the solid red line to the dotted red line in Fig 10C). The greater effect of cost on mixtures means the advantage of sequences over mixtures in Fig 10A (and Fig 6A) is removed in Fig 10C. With costs set to 0 time-to-resistance for mixture divided by sequence is 0.8 in Fig 10A and 1 in Fig 10C. At higher insecticide effectiveness the benefit of mixture over sequence is improved as seen by comparing Fig 10D to 10B.

A summary of the effects of inputs on the evolution of resistance for mixtures is given in Table 3 and for the difference between mixtures and sequences in Table 4. The mechanisms identified in the tables are covered further in the discussion.

Table 3. Effect of inputs on resistance when insecticides used in a mixture

Parameter to increase | effect on resistance evolution | Mechanism ----------------------- | ---------- | --------------------------------------------- 1. Effectiveness | slower | one insecticide reduces the fitness of individuals resistant to the other thereby reducing selection pressure for the other 2. Exposure | faster (but less than for insecticide used alone) | reduced fitness of individuals susceptible to one insecticide increases selection pressure for that resistance. However at the same time selection pressure is reduced by reduced fitness of resistant individuals caused by the other insecticide 3. Dominance of resistance | faster | increased fitness of heterozygotes 4. Resistance restoration | faster | increased fitness of resistants 5. Frequency | faster | less change needed to reach resistance threshold 6. Cost of resistance | slower | reduced fitness of resistants in absence of insecticide

Table 4. Effect of inputs on the difference between mixture and sequential use

Parameter to increase | increase favours mix or sequence | Mechanism ----------------------- | -------- | --------------------------------------------- 1. Effectiveness | mixture | Higher effectiveness gives faster resistance for sequence and slower resistance in mixture 2. Exposure | sequence | Higher exposure gives faster resistance for sequence and mixture but the greater effect on mixture favours sequence. 3. Dominance of resistance | neither | Higher dominance gives faster resistance in both sequences and mixtures such that the difference between them is not changed. 4. Resistance restoration | neither | As for dominance of resistance. 5. Frequency | neither | As for dominance and resistance restoration. 6. Cost of resistance | mixture | Higher costs slow the evolution of resistance more in a mixture than when used in sequence (but with higher costs in a sequence there is a greater chance for resistance levels to decline for the insecticide not being used).

Linkage disequilibrium

Positive linkage disequilibrium is said to exist when the frequency of two alleles at different loci is higher than expected if they were just selected individually. It has previously been suggested [@Curtis1985][@WHO2012] that using an insecticide to which there was existing resistance in a mixture with a new insecticide with little resistance could lead to more rapid evolution of resistance to the latter through linkage disequilibrium. We explicitly look at the example of a mixture of a new insecticide and one with pre-existing resistance (Figures 9 C&D). In both cases resistance to the newer insecticide with lower starting resistance (in red) increases faster when it is used alone (red dotted) than when it is used in a mixture (red solid). This suggests that linkage disequilibrium does not restrict the potential benefits of mixtures in all situations. These results support the statement in GPIRM [@WHO2012] that : "Even though the risk for linkage disequilibrium exists, mixtures may still be the most attractive IRM tool …". The dynamics of LD is too complex to discuss in detail here but interested readers can find more discussion in [@Curtis1985][@WHO2012][@Levick2017].

Setting model inputs from field data

The main model inputs quantifying properties of the mosquitoes and insecticides can be derived from data on the relative survival of the different genotypes (RR,SR,SS) as described in Table 5.

Table 5. Equations to calculate model input parameters from field data. Where RRfit, SRfit, SSfit are fitnesses of genotypes. Derived from [@Levick2017].

Parameter | Calculation ---------------------- | ---------------------------------------------------- 1. Effectiveness | 1 - SSfit in presence of insecticide 2. Exposure | estimated 3. Resistance restoration | (RRfit-SSfit) / Effectiveness in presence of insecticide 4. Dominance of resistance | (SRfit-SSfit)/(RRfit-SSfit) in presence of insecticide 5. Frequency | estimated 6. Cost of resistance | 1 - RRfit in absence of insecticide 7. Dominance of cost | (SRfit-RRfit)/(RRfit-SSfit) in absence of insecticide

Exposure, the proportion of mosquitoes that are exposed to the insecticide, is a property of the location and will depend on the use of nets or IRS and mosquito behaviour. The frequency of resistance alleles can be measured, although when at low levels a large number of mosquitoes would need to be sampled to detect low frequencies.

As an example of making these calculations we use data for Anopheles gambiae on pyrethroid (Alpha-cypermethrin) resistance associated with the Kdr mutation [@Kolaczinski2000] and carbamate (Bendiocarb) resistance associated with the Ace1 mutation [@Essandoh2013]. This is not intended to be a definitive prediction of resistance evolution, sample sizes for some genotypes are low so fitnesses are approximate. It is intended to demonstrate how models of resistance evolution can be set to field data and hence bring the two approaches closer together.

Data were provided as numbers alive and dead by genotype in both publications. We used these to calculate genotype-specific survival values as an indication of fitness. We used the non-exposed control data for the pyrethroid [@Kolaczinski2000] to estimate the cost of resistance and to rescale the fitness estimates on a scale of 0-1 where 1 is the fitness of the unexposed SS (i.e. we divided all pyrethroid survivals by that of the unexposed SS genotype, 0.83). In the absence of unexposed data for the carbamate [@Essandoh2013] we assumed no costs of resistance and that survival of the non-exposed SS would be 1. From these estimates we calculated the values of our model inputs (Table 6).

Table 6. Model input parameters calculated from field data.

Field data|Pyrethroid / Kdr[@Kolaczinski2000]|rescaled Pyrethroid / Kdr[@Kolaczinski2000]|Carbamate / Ace1[@Essandoh2013] ------------------------------|-------|-------|------- exposed survival RR|0.83|1|0.84 exposed survival SR|0.72|0.87|0.56 exposed survival SS|0.59|0.71|0.02 unexposed survival RR|0.67|0.81||- unexposed survival SR|0.50|0.60|- unexposed survival SS|0.83|1|- Calculated model inputs| | | Effectiveness|0.41|0.29|0.98 Resistance restoration|0.60|1|0.84 Dominance of resistance|0.52|0.54|0.66 Cost of resistance|-|0.19|- Dominance of cost|-|1*|-

*For the pyrethroid the unexposed survival of the SR (0.5) was less than that for the RR (0.67) This indicates underdominance which is biologically unlikely for unexposed genotypes so we set dominance at 1 so that SR and RR genotypes have the same fitness costs. The authors indicated that suprisingly high mortality in the control huts may have been due to rough handling [@Kolaczinski2000] and this could contribute to the unexpected values here.

The results of using these calculated inputs in one model example are shown in Figure 11. Remember that this is a demonstration exercise for setting model parameters from field data and not a definitive prediction for whether using a pyrethroid and a carbamate in a mixture or sequence is likely to be preferable.

Discussion

Resistance responds to insecticide use in the model in ways that we can explain mechanistically through the process of selection. This mechanistic explanation can help us develop a more robust understanding of how we would expect resistance to evolve in the field. We start by explaining the evolution of resistance to a single insecticide, which is relatively straightforward, and move on to explaining the effect of two insecticides in sequences and mixtures.

When single insecticide use was investigated (Figs 4-5) higher values of parameters that increased the selective advantage of resistance all led to more rapid evolution of resistance, whereas higher values for cost led to slower evolution of resistance (Fig 5B). Mechanisms explaining the evolution of resistance to single insecticide use, are summarised in Table 2. Referring to Figure 2 can help explain these responses. Exposure and effectiveness act by decreasing the fitness of SS and SR more than RR, dominance of resistance acts by increasing fitness of SR and resistance-restoration acts by increasing fitness of RR & SR. Cost acts by decreasing the fitness of RR and SR among those mosquitoes not exposed.

Two insecticides

Responses to input parameters are different when two insecticides are used either in sequence or in mixtures as compared to when used alone. These differences can help us to understand the dynamics of resistance evolution. In each panel of Figures 6-10 the solid lines, indicating that an insecticide is being used at the same time as another in a mixture, are always shallower than the same colour dashed line indicating the insecticide is being used in isolation as part of a sequence. The shallower curves for the mixture can be explained by each insecticide killing individuals that are resistant to the other insecticide. Because resistant individuals are being killed, the selection pressure for that resistance is lower. Thus each insecticide reduces the selection pressure for resistance to the other and could be said to 'protect' the other. This protection means that evolution of resistance to one insecticide is always slower when it is used in a mixture than when used alone. This protection has been termed multiple intra-generational killing [@Consortium2013] because individuals have the potential to be killed by more than one insecticide.

Whilst this protection ensures that resistance to each insecticide evolves more slowly in a mixture it does not guarantee that a mixture strategy is preferable to a sequence once the resistance to both insecticides is taken into account. In a mixture, selection pressure is applied by both insecticides at the same time, whereas in a sequence there is no selection for the insecticide not being used. Resistance can evolve more slowly for both insecticides in sequential use, despite being faster for each, because they occur one after the other. In Figure 6A it takes 60 generations for resistance to reach the threshold for each insecticide in the sequence, and 90 when they are applied together in the mixture, so the sequence is still favoured because 2x60 is greater than 90.

This advantage of sequence over mixtures can be removed by increasing the effectiveness of one of the insecticides. Increasing the effectiveness of one insecticide (Fig 6B) decreases the total time for the sequence from 120 to 90 and increases that for the mixture from 90 to 115 thus changing the order and favouring the mixture. The mechanism for this change can be seen by comparing the shape of the curves in Figs 6A & B. In the mixture, resistance to insecticide 1 (with the increased effectiveness) rises at a similar speed in both figures. However, the increased effectiveness of insecticide 1 causes resistance to insecticide 2 to increase more slowly initially in the mixture. Once insecticide 1 reaches resistance of around 50% at around 50 generations, the curve for insecticide 2 becomes steeper. Once the first insecticide reaches the resistance threshold of 50% it kills fewer individuals that are resistant to the second insecticide, it loses its protective effect and thus stops slowing the rise in resistance to the second insecticide. A similar effect is visible when the effectiveness of both insecticides are increased (Fig 7C). The identical resistance curves for each insecticide in the mixture are shallower than at the lower effectiveness (Fig 6A) because more individuals resistant to each insecticide are being killed by the other insecticide. In this case the 'protection' given to both insecticides by the other declines at the same rate so we do not see the change in slope seen in Fig 6B. In both cases increasing the effectiveness of insecticides increases their protective effect when used in mixtures and favours mixtures over sequences.

Fig 6B also demonstrates that when two insecticides have a different effectiveness in a mixture, other parameters being equal, resistance will be expected to increase faster to the more effective insecticide (the red curve in this case). The more effective insecticide prompts both a) higher selection pressure for its own resistance and b) greater 'protection' reducing the rise in resistance to the other insecticide. The more effective insecticide has a faster evolution of resistance than the less effective one but this is not caused by the presence of the less effective one. Indeed the presence of the less effective one is still slowing the evolution of resistance to the more effective one, as indicated by the red solid curve for use in the mixture being shallower than the red dashed curve for use alone as part of the sequence.

Difference between exposure and effectiveness

It seems initially counter-intuitive that whereas increasing effectiveness slows the development of resistance for mixtures, increasing exposure speeds it up. Increasing the exposure to both insecticides results in a decrease in time-to-resistance (from 90 generations in Fig 6A to 40 generations in Fig 6C) where increasing the effectiveness of both insecticides by the same amount results in a slight increase (from 90 in Fig 6A to 115 in Fig 6B). This contrasts with the identical effects that exposure and effectiveness have on a single insecticide (Figs 4 A,B) and on insecticides in sequence (the red dashed lines are the same in Figs 6B & C). This difference, that exposure reduces times-to-resistance for a mixture and effectiveness increases it, was unexpected for us and we are keen to explore it in more detail. The most likely explanation is that increasing exposure increases overall selection pressure and lowers time to resistance, whereas increasing effectiveness increases mutual protection of the insecticides in a mixture and lowers the speed at which resistance evolves.

Costs of resistance

Fitness costs of resistance have been widely reported(review in [@Kliot2012]) but it has been suggested that few of these are directly relevant to control failure in the field [@Ffrench-Constant2017]. Because of this likely rarity of costs we did not include them in most of our analyses (with the exception of Fig 10).

Costs of resistance favour mixtures in these analyses because the evolution of resistance is slowed more by cost in a mixture than in a sequence (compare the red solid lines for a mixture to the red dotted line for a sequence in Fig 10C and D). It makes sense that the protective effect of the other insecticide in a mixture combines with the cost of resistance to reduce selection pressure more than in a sequence.

The resistance curves for two insecticides in a mixture that differ only in their cost of resistance, diverge more as the frequency of resistance increases (Fig 10C). This makes intuitive sense, as the frequency of resistance genes increases the impact of fitness costs associated with those genes also increases. Thus differences due to cost are less noticeable at low resistance frequencies and become more noticeable at higher resistance frequencies. Thus if a lower resistance threshold was used cost would have less of an effect and if a higher resistance threshold was used cost would have more of an effect. Also note that if the costs are set even higher than we set them here (e.g. above 0.2) they can maintain resistance frequencies below 1 or lead to their decline. However, we do not expect such high costs to occur in operational settings.

Declines in resistance frequencies are likely to occur when there are costs of resistance and the insecticide is not being used [@Ffrench-Constant2017]. Thus the frequency of resistance for the first insecticide in a sequence, after it has been replaced, may decline to a level where the first insecticide becomes operationally useful again. Consequently, a repeated sequential strategy could offer benefits over mixtures not considered here. To quantify these benefits one would need to use a more detailed metric than the time-to-resistance we have used. We quantified the number of generations until resistance thresholds are first reached. Other measures for quantifying resistance over time can be imagined, for example the mean resistance or proportion of generations below a resistance threshold.

Using field derived inputs

Model behaviour when using the field derived inputs (Fig. 11) can be explained using the understanding developed from the preceeding runs. Resistance to the carbamate (blue) rises quickly partly because it has high effectiveness (0.98). Resistance to the pyrethroid (red) increases more slowly due to it's much lower effectiveness (0.29) and high cost of resistance. In the mixture the resistance to the pyrethroid declines at the start when resistance to the carbamate is low. This illustrates the protective effect of the insecticide with the higher effectiveness. That protective effect causes the mixture to out-perform the sequence in this example. Indeed, if exposure is set lower to 0.5 then resistance to the pyrethroid declines further and the resistance threshold is not reached. This is a combination of the low effectiveness creating a low selection pressure for the resistance and the high cost creating a selection pressure against it. The low effectiveness and high cost are partly due to the surprisingly high mortality in unexposed mosquitoes in the field data we used [@Kolaczinski2000].

Consistency with previous work

The predicted responses of resistance presented here can be put into the context of previous work and recommendations. A recent review [@Consortium2013] found that 8 out of 10 empirical studies and 11 out of 14 modelling studies favoured mixtures over sequences. In the 3 other modelling studies the relative performance of mixtures and sequences were dependent on the values of other inputs as we also have shown.

The review concluded that the advantage of mixtures is greatest when seven other conditions are satisfied [@Consortium2013] : i) rare starting resistance , ii) independent loci (i.e. no cross resistance), iii) high recombination, iv) high susceptible mortality, v) resistance is recessive, vi) similar persistence of insecticides, and vii) some of population remains untreated. Their points iv) and vii) agree with our results that mixtures are favoured by high effectiveness and low exposure. Some of the other conditions, notably i) the starting frequency of resistance and v) the dominance of resistance did not alter whether mixtures or sequences were favoured in our model [@Levick2017].

Earlier modelling work has made sometimes contrasting predictions and these can be explained by the details of their structure and input choices. One model predicted that "the use of mixtures is always more effective in delaying the onset of resistance often by many orders of magnitude"[@Mani1985]. That model assumed that SS genotypes were always killed, RR always survived, and that a constant 10% of mosquitoes escaped the insecticide. This would be represented in our model by a very high 'insecticide effectiveness' of 1, a 'resistance restoration' of 1 and an 'exposure' of 0.9. This result is consistent with our prediction that the very high effectiveness is likely to favour mixtures over sequences (e.g. see Fig 7D and Table 4).

A subsequent model lead to the following more pessimistic conclusion about the potential of mixtures : “As a result of incomplete coverage and residue decay, the mortality of susceptible homozygotes is rarely consistently high enough for pesticide mixtures to be effective” [@Roush1998]. This "mortality of susceptible homozygotes" is equivalent to our insecticide effectiveness. Again this points to the importance of effectiveness and the question of whether the effectiveness levels needed to make mixtures better are attainable operationally. Consistent with our results [@Roush1998] also showed (Figure 3) that decreasing exposure favours mixtures over sequences.

Caveats

There is evidence of the involvement of multiple genes in resistance mechanisms[@Ffrench-Constant2004][@Ffrench-Constant2013][@Gressel2017]. Our model, in common with almost all previous ones [@Bourguet2010], assumes a single ‘resistance’ gene per insecticide. In many cases we believe this is reasonable. Target-site resistance has been shown to be mostly determined by single genes, and control failures in general have mostly been attributed to the effect of single major genes [@Roush1989][@Ffrench-Constant2004][@Ffrench-Constant2013]. One view is that control doses of insecticides in the field are likely to lead to the selection of single major genes whereas the lower doses used to select laboratory strains are more likely to lead to polygenic resistance [@Ffrench-Constant2013]. There is also evidence that low insecticide doses in agriculture have led to polygenic resistance [@Gressel2017]. To represent polygenic resistance a different 'quantitative genetic' modelling approach would be required [@Via1986][@Firko1990]. A recent review of 187 models of the evolution of resistance [@Bourguet2010] did not find any that represented quantitative multiple gene resistance. We are only aware of limited published work applying a quantitative genetic approach to assessing insecticide use strategies in the context of polygenic resistance [@Via1986][@Gardner1998][@Gardner1999]. If it becomes clearer that polygenic resistance is an issue for vector control then there will be a need for appropriate quantitative genetic models. We would be keen to contribute to that development.

Poor implementation

One aspect of the results presented here can seem to contradict published advice and has been questioned when we have presented model results. Our model predicts that lower insecticide effectiveness or exposure will lead to slower spread of resistance in sole use and sequences (Figs 4 & 6, Table 2). This is consistent with expectations that lower effectiveness and exposure would lead to lower mortality of susceptibles and thus lower selection pressure for resistance [e.g. @Roush1998]. This can seem to contradict recommendations that poor implementation of insecticide interventions can promote the development of resistance. e.g. from [@FAO2012] "Certain pest control practices have consistently been shown to exacerbate the loss of susceptible pest populations and the development of resistance. These include: ... - the use of application rates that are below or above those recommended on the label, - poor coverage of the area being treated ...".

Poor implementation would likely cause lower effectiveness of, and exposure to, insecticides, so how can this increase the potential for resistance when our model suggests it would be decreased ? There are two main potential explanations for this apparent difference, a) dominance and b) polygenic resistance. The effect of dominance is that poor application can reduce the mortality of heterozygotes (SR) thus effectively increasing the dominance of the resistant allele and increasing its rate of evolution [@Bourguet2000][@Barbosa2011][@FAO2012] see Fig 1 of [@Levick2017]. Our model results are consistent with this, showing that increasing dominance of resistance leads to faster spread of resistance (Fig 4C, 8). So, the effect of poor implementation of an insecticide intervention on resistance could be either positive or negative depending on the relative effects on dominance of the resistance gene, insecticide effectiveness and exposure. Similarly, the effects of poor application on the relative benefits of mixtures and sequences will depend on these trade-offs because of our result that higher effectiveness favours mixtures, higher exposure favours sequences and higher dominance favours neither (Table 4). More research on how the mortality of different genotypes (susceptible, resistant and heterozygous) is affected by poor application rates or insecticide degradation over time would help make this clearer. Secondly, if resistance is predominantly polygenic rather than monogenic poor implementation is likely to promote resistance by allowing the build up of many genes of small effect [@Gressel2017]. Indeed modelling work from 15 years ago [@Gardner1998][@Gardner1999] proposes an alternation of low and high doses with the low doses limiting the evolution of monogenic resistance and the high doses limiting that of the polygenic resistance. More research on the relative importance of monogenic and polygenic resistance for vector control would help here.

Experience from agriculture where single interventions have been associated with a rapid development of resistance have led for calls for a more Integrated Vector Management composed of a series of partially effective tools [@Thomas2012]. It is suggested that such an integrated approach is more sustainable and 'evolution-proof' [@Thomas2012]. A similar combining of interventions has been advocated to cover mosquitoes exhibiting different behaviours [@Killeen2017]. Mosquitoes with different behaviours e.g. feeding outdoors and/or on animals are likely to result in lower exposures to insecticides in nets or sprayed on walls. This lower exposure could favour mixtures over sequences as we have indicated. In both cases our modelling approach could be used to investigate implications for the development of resistance.

Even if the assumption that single genes are responsible for control failures is correct there is still considerable other uncertainty around the development of insecticide resistance in operational field settings [@Ranson2016]. Our model is able to represent much of that uncertainty, yet in this paper we restrict our explorations to a limited region of parameter space. Here we hold most inputs constant at intermediate values while modifying single inputs in isolation. This is a deliberate tactic to develop understanding and communicate the key mechanisms. There is the potential that the system may behave differently in the field. However, the results presented here are supported by an earlier analysis where we ran 10,000 unique scenarios varying more inputs [@Levick2017]. This previous analysis also highlighted the relative importance of insecticide effectiveness and exposure in determining the relative performance of mixtures and sequences.

In addition we provide an on-line user interface to the model (Fig 3) [@South2017a] enabling readers with no coding experience to investigate model behaviour in areas of parameter space we have not shown. The user-interface allows inputs to be set for two scenarios and results to be viewed side by side. For example you can set a higher insecticide effectivenesses in one scenario to see how this favours mixtures over sequences (Fig 3) or increase the exposure levels to see how this favours sequences.

There are other model behaviours that can be changed by modifying model inputs within the computer code. For example, here we assume that insecticide interaction in a mixture is multiplicative, e.g. if the probability of surviving exposure to insecticide A alone is 0.3 and of surviving insecticide B alone is 0.2, then the probability of surviving exposure to a mixture of A and B is 0.3*0.2 = 0.06. This is consistent with recent work on aphids [@Taillebois2016], although that did also indicate some synergistic effects. Modifying model inputs would allow us to incorporate such synergistic or antagonistic effects of mixtures, i.e. that insecticides used together perform better or worse than when used separately. This would allow the model to represent cross-resistance where a resistance gene provides resistance to more than one insecticide. Similarly behavioural resistance is likely to offer resistance against more than one insecticide [@WHO2012] but until there are more data on whether behavioural resistance is genetic and heritable we would not recommend representing it in this way.

In addition to these uncertainties regarding resistance it is important to remember that our model examines solely the evolution of resistance and not the control of mosquitoes or malaria. Different insecticide use strategies, as well as affecting the evolution of resistance, will also influence both the numbers of vectors and their propensity to transmit the disease [@Rivero2010]. There are also potential knock-on effects of resistance itself; reduced vector transmission of disease (e.g. by reduced lifespan of the vector) or increased transmission (e.g. by increasing the susceptibility of the vector to the disease) [@Rivero2010]. This model can explore when a mixture strategy may be favoured over a sequence purely for the evolution of resistance. Other operational factors also need to be considered. On the positive side a mixture is likely to kill more vectors and may also limit their ability to transmit the disease. On the negative side a mixture may be more expensive. Other attributes of insecticides such as their repellency are also relevant. Insecticides with some repellency will reduce exposure to the insecticide and thus, as shown here, are likely to slow the evolution of resistance and favour mixtures over sequences. However the reduced exposure resulting from repellency also reduces the population protection provided by the insecticide [@Birget2015].

Our model would predict that the best strategy for limiting the development of resistance in the short term would be to use no insecticide. Of course, using no insecticide is likely to have serious consequences for mosquito and malaria abundance. The challenge is to develop good strategies for delivering mosquito and malaria control in the short term while sustaining their operational capacity in the longer term by restricting the development of resistance.

Conclusions

We have developed a flexible, modern model of the evolution of insecticide resistance that can help us understand the mechanisms by which different insecticide application strategies may be favoured. It replicates previous modelling work and has the potential to investigate further operational strategies for established and new insecticides. For assessing the relative benefits of mixtures and sequences we have demonstrated that mixtures are likely to be favoured by high insecticide effectiveness, low exposure and high costs of resistance. Similar processes will occur in other strategies such as rotations and mosaics. This accessible modelling approach can help to promote discussion and exploration of the likely effect on resistance evolution of different insecticide use strategies in tandem with the collection of empirical data to refine and test predictions.

Declarations

Ethics approval and consent to participate

Not applicable

Consent for publication

Not applicable

Availability of data and material

Project name: resistance Project home page: https://github.com/AndySouth/resistance Archived version: 10.5281/zenodo.889012 Operating system: Platform independent Programming language: R License: GNU General Public License v3.0

Competing interests

AS received a small grant from Bayer S.A.S. to support a few months work and attendance of the American Society of Tropical Medicine and Hygiene Annual Meeting.

Funding

This work was supported by: Bill and Melinda Gates Foundation, grant #37999.01 to IMH; and Wellcome Trust Institutional Strategic Strengthening fund grant #105620/Z/14/Z to Liverpool School of Tropical Medicine; sub-funding to IMH; a small grant from Bayer S.A.S. to AS; and IVCC funding to IMH.

Authors' contributions

AS conceived the paper, re-factored the model code, performed the analyses and wrote the first and subsequent versions of the paper. IH conceived the modelling approach, developed the first version of the model in collaboration with Bethany Levick, explained concepts to AS and reviewed and commented on all manuscript versions.

Acknowledgements

Thank you to Bethany Levick who developed the first version of the model code in collaboration with IH. Thank you to David Malone, Justin McBeath Hilary Ranson, Sarah Rees, Graham Small and an anonymous reviewer for commenting constructively on an earlier draft of this paper. Thank you to Brian Bourke and phylopic.org for the mosquito image used in Figure 1.

\pagebreak

Figures

plot_fit_calc( effectiveness=0.8, resistance_restoration=0.5, dominance_restoration=0.7, dominance_cost=0.7, cost=0.3 )

  library(cowplot)

  #list to put plots in for cowplot
  plotlist <- list(4)

  #A effectiveness  
  a <- setExposure(exposure=0.5, insecticideUsed='insecticide1')

  i1 <- setInputOneScenario( max_gen = 500,
                      h.RS1_A0 = 0.5, 
                      h.RS2_0B = 0.5,
                      a = a,
                      phi.SS1_A0 = 0.2,
                      phi.SS2_0B = 0.5,
                      rr_restoration_ins1 = 0.5,
                      rr_restoration_ins2 = 0.5 )

  i2 <- setInputOneScenario( max_gen = 500,
                      h.RS1_A0 = 0.5, 
                      h.RS2_0B = 0.5,
                      a = a,
                      phi.SS1_A0 = 0.5,
                      phi.SS2_0B = 0.5,
                      rr_restoration_ins1 = 0.5,
                      rr_restoration_ins2 = 0.5 )

  i3 <- setInputOneScenario( max_gen = 500,
                      h.RS1_A0 = 0.5, 
                      h.RS2_0B = 0.5,
                      a = a,
                      phi.SS1_A0 = 0.8,
                      phi.SS2_0B = 0.5,
                      rr_restoration_ins1 = 0.5,
                      rr_restoration_ins2 = 0.5 )

  #input <- cbind(input, inputOneScenario)
  input <- cbind(i1,i2,i3)

  listOut <- runModel2( input )

  df_resist <- get_resistance(locus=1, listOut)

  #A  
  plotlist[[1]] <- ggplot(df_resist, aes(x=generation, y=resistance, colour=factor(phi.SS1_A0))) + 
         theme_bw() +  
         theme(legend.position = "bottom", legend.key = element_blank()) + 
         guides(colour = guide_legend(reverse=TRUE)) +  
         labs(colour = "effectiveness") +
         coord_trans(y = "log10") +    
         geom_line()   


  #B exposure
  exposure <- 0.2
  a <- setExposure(exposure = exposure, insecticideUsed='insecticide1')

  i1 <- setInputOneScenario( max_gen = 500,
                      h.RS1_A0 = 0.5, 
                      h.RS2_0B = 0.5,
                      a = a,
                      exposure = exposure, #used for post-run analysis
                      phi.SS1_A0 = 0.5,
                      phi.SS2_0B = 0.5,
                      rr_restoration_ins1 = 0.5,
                      rr_restoration_ins2 = 0.5 )

  exposure <- 0.5
  a <- setExposure(exposure = exposure, insecticideUsed='insecticide1')   

  i2 <- setInputOneScenario( max_gen = 500,
                      h.RS1_A0 = 0.5, 
                      h.RS2_0B = 0.5,
                      a = a,
                      exposure = exposure, #used for post-run analysis
                      phi.SS1_A0 = 0.5,
                      phi.SS2_0B = 0.5,
                      rr_restoration_ins1 = 0.5,
                      rr_restoration_ins2 = 0.5 )

  exposure <- 0.8
  a <- setExposure(exposure = exposure, insecticideUsed='insecticide1')  

  i3 <- setInputOneScenario( max_gen = 500,
                      h.RS1_A0 = 0.5, 
                      h.RS2_0B = 0.5,
                      a = a,
                      exposure = exposure, #used for post-run analysis
                      phi.SS1_A0 = 0.5,
                      phi.SS2_0B = 0.5,
                      rr_restoration_ins1 = 0.5,
                      rr_restoration_ins2 = 0.5 )

  input <- cbind(i1,i2,i3)

  listOut <- runModel2( input )

  df_resist <- get_resistance(locus=1, listOut)

  #B  
  plotlist[[2]] <- ggplot(df_resist, aes(x=generation, y=resistance, colour=factor(exposure))) + 
         theme_bw() +
         theme(legend.position = "bottom", legend.key = element_blank()) + 
         guides(colour = guide_legend(reverse=TRUE)) +
         labs(colour = "exposure", y=NULL) +
         coord_trans(y = "log10") +   
         geom_line()   

  #plot_grid( plotlist[[1]],plotlist[[2]],ncol=2, rel_widths=c(1,1), labels='AUTO') 


  #C single insecticide dominance
  a <- setExposure(exposure=0.5, insecticideUsed='insecticide1')

  i1 <- setInputOneScenario( max_gen = 500,
                      h.RS1_A0 = 0.2, 
                      h.RS2_0B = 0.5,
                      a = a,
                      phi.SS1_A0 = 0.5,
                      phi.SS2_0B = 0.5,
                      rr_restoration_ins1 = 0.5,
                      rr_restoration_ins2 = 0.5 )

  i2 <- setInputOneScenario( max_gen = 500,
                      h.RS1_A0 = 0.5, 
                      h.RS2_0B = 0.5,
                      a = a,
                      phi.SS1_A0 = 0.5,
                      phi.SS2_0B = 0.5,
                      rr_restoration_ins1 = 0.5,
                      rr_restoration_ins2 = 0.5 )

  i3 <- setInputOneScenario( max_gen = 500,
                      h.RS1_A0 = 0.8, 
                      h.RS2_0B = 0.5,
                      a = a,
                      phi.SS1_A0 = 0.5,
                      phi.SS2_0B = 0.5,
                      rr_restoration_ins1 = 0.5,
                      rr_restoration_ins2 = 0.5 )

  input <- cbind(i1,i2,i3)

  listOut <- runModel2( input )

  df_resist <- get_resistance(locus=1, listOut)

  #C 
  plotlist[[3]] <-ggplot(df_resist, aes(x=generation, y=resistance, colour=factor(h.RS1_A0))) + 
         theme_bw() +
         theme(legend.position = "bottom", legend.key = element_blank()) + 
         guides(colour = guide_legend(reverse=TRUE)) + 
         labs(colour = "dominance of resistance") +
         coord_trans(y = "log10") + 
         geom_line()   


  #D single insecticide resistance restoration 

  a <- setExposure(exposure=0.5, insecticideUsed='insecticide1')

  i1 <- setInputOneScenario( max_gen = 500,
                      h.RS1_A0 = 0.5, 
                      h.RS2_0B = 0.5,
                      a = a,
                      phi.SS1_A0 = 0.5,
                      phi.SS2_0B = 0.5,
                      rr_restoration_ins1 = 0.2,
                      rr_restoration_ins2 = 0.5 )

  i2 <- setInputOneScenario( max_gen = 500,
                      h.RS1_A0 = 0.5, 
                      h.RS2_0B = 0.5,
                      a = a,
                      phi.SS1_A0 = 0.5,
                      phi.SS2_0B = 0.5,
                      rr_restoration_ins1 = 0.5,
                      rr_restoration_ins2 = 0.5 )

  i3 <- setInputOneScenario( max_gen = 500,
                      h.RS1_A0 = 0.5, 
                      h.RS2_0B = 0.5,
                      a = a,
                      phi.SS1_A0 = 0.5,
                      phi.SS2_0B = 0.5,
                      rr_restoration_ins1 = 0.8,
                      rr_restoration_ins2 = 0.5 )

  input <- cbind(i1,i2,i3)

  listOut <- runModel2( input )

  df_resist <- get_resistance(locus=1, listOut)

  #D  
  plotlist[[4]] <- ggplot(df_resist, aes(x=generation, y=resistance, colour=factor(rr_restoration_ins1))) + 
         theme_bw() +
         theme(legend.position = "bottom", legend.key = element_blank()) + 
         guides(colour = guide_legend(reverse=TRUE)) +
         labs(colour = "resistance restoration", y=NULL) +
         coord_trans(y = "log10") +   
         geom_line()   

  cowplot::plot_grid( plotlist[[1]],plotlist[[2]],plotlist[[3]],plotlist[[4]], ncol=2, rel_widths=c(1,1), labels='AUTO')

  #single insecticide frequency fig

  #A starting frequency
  a <- setExposure(exposure=0.5, insecticideUsed='insecticide1')

  i1 <- setInputOneScenario( max_gen = 500,
                      P_1 = 0.001,       
                      h.RS1_A0 = 0.5, 
                      h.RS2_0B = 0.5,
                      a = a,
                      phi.SS1_A0 = 0.5,
                      phi.SS2_0B = 0.5,
                      rr_restoration_ins1 = 0.5,
                      rr_restoration_ins2 = 0.5 )

  i2 <- setInputOneScenario( max_gen = 500,
                      P_1 = 0.01,       
                      h.RS1_A0 = 0.5, 
                      h.RS2_0B = 0.5,
                      a = a,
                      phi.SS1_A0 = 0.5,
                      phi.SS2_0B = 0.5,
                      rr_restoration_ins1 = 0.5,
                      rr_restoration_ins2 = 0.5 )

  i3 <- setInputOneScenario( max_gen = 500,
                      P_1 = 0.1,       
                      h.RS1_A0 = 0.5, 
                      h.RS2_0B = 0.5,
                      a = a,
                      phi.SS1_A0 = 0.5,
                      phi.SS2_0B = 0.5,
                      rr_restoration_ins1 = 0.5,
                      rr_restoration_ins2 = 0.5 )

  input <- cbind(i1,i2,i3)

  listOut <- runModel2( input )

  df_resist <- get_resistance(locus=1, listOut)

  #A  
  plotlist[[1]] <-  ggplot(df_resist, aes(x=generation, y=resistance, colour=factor(P_1))) + 
         theme_bw() +
         theme(legend.position = "bottom", legend.key = element_blank()) + 
         guides(colour = guide_legend(reverse=TRUE)) + 
         labs(colour = "start\nfrequency") +
         coord_trans(y = "log10") + 
         geom_line()   

  #B changing costs at dominance 0.5

  a <- setExposure(exposure=0.5, insecticideUsed='insecticide1')

  i1 <- setInputOneScenario( max_gen = 500,
                      P_1 = 0.01,       
                      h.RS1_A0 = 0.5, 
                      h.RS2_0B = 0.5,
                      a = a,
                      phi.SS1_A0 = 0.5,
                      phi.SS2_0B = 0.5,
                      rr_restoration_ins1 = 0.5,
                      rr_restoration_ins2 = 0.5,
                      z.RR1_00 = 0,
                      h.RS1_00 = 0.5)

  i2 <- setInputOneScenario( max_gen = 500,
                      P_1 = 0.01,       
                      h.RS1_A0 = 0.5, 
                      h.RS2_0B = 0.5,
                      a = a,
                      phi.SS1_A0 = 0.5,
                      phi.SS2_0B = 0.5,
                      rr_restoration_ins1 = 0.5,
                      rr_restoration_ins2 = 0.5,
                      z.RR1_00 = 0.1,
                      h.RS1_00 = 0.5 )

  i3 <- setInputOneScenario( max_gen = 500,
                      P_1 = 0.01,       
                      h.RS1_A0 = 0.5, 
                      h.RS2_0B = 0.5,
                      a = a,
                      phi.SS1_A0 = 0.5,
                      phi.SS2_0B = 0.5,
                      rr_restoration_ins1 = 0.5,
                      rr_restoration_ins2 = 0.5,
                      z.RR1_00 = 0.2,
                      h.RS1_00 = 0.5 )

  input <- cbind(i1,i2,i3)

  listOut <- runModel2( input )

  df_resist <- get_resistance(locus=1, listOut)

  #B  
  plotlist[[2]] <- ggplot(df_resist, aes(x=generation, y=resistance, colour=factor(z.RR1_00))) + 
         theme_bw() +
         theme(legend.position = "bottom", legend.key = element_blank()) + 
         guides(colour = guide_legend(reverse=TRUE)) + 
         labs(colour = "cost of\nresistance") +
         coord_trans(y = "log10") + 
         geom_line()   

  cowplot::plot_grid( plotlist[[1]],plotlist[[2]], ncol=2, rel_widths=c(1,1), labels='AUTO')  
  #cowplot::plot_grid( plotlist[[1]],plotlist[[2]],plotlist[[3]],plotlist[[4]], ncol=2, rel_widths=c(1,1), labels='AUTO')

#4 panel figure looking at effectiveness and exposure


#test of rearranging figures on a grid
#eff1  0.5  0.8
#exp
#0.5    1   2
#0.8    3   4

#layout(matrix(1:4,2,2, byrow = TRUE))

#add lower panel for legend
#layout(matrix(c(1,2,3,4,5,5),3,2, byrow = TRUE), heights=c(1,1,0.5), respect=FALSE)

#layout( matrix(c(1,2,3,4,5,5),3,2, byrow = TRUE), heights=c(1,1,0.2), 
# respect=matrix(c(1,1,1,1,0,0),3,2, byrow = TRUE) )
#layout.show(n=5)

# add side panels for labels
layout( matrix(c(1:6,0,7,7,0,8,9),4,3, byrow = TRUE), heights=c(1,1,0.15,0.2), widths=c(0.1,1,1) )
#layout.show(n=9)

#effectiveness on mixture


#label
par(mar=c(0, 0, 0, 0), pty='m') #b,l,t,r default c(5, 4, 4, 2)
plot( 0, type="n", axes=FALSE, ann=FALSE, pty='m')
mtext('exposure = 0.5', line= -2, side=2, cex=0.9) #side=1b,2l,3t,4r

#base scenario

par(mar=c(1.5, 4, 2, 0)) #b,l,t,r default c(5, 4, 4, 2)

runcurtis_f2( max_gen=500,  P_1 = 0.01 , P_2 = 0.01 , h.RS1_A0 = 0.5 , h.RS2_0B = 0.5 , exposure = 0.5 , phi.SS1_A0 = 0.5 , phi.SS2_0B = 0.5 , rr_restoration_ins1 = 0.5 , rr_restoration_ins2 = 0.5 , addCombinedStrategy = FALSE, strategyLabels = c('s','','a','m'), xlabs = FALSE, cex.axis = 0.8, addLegend=FALSE, main='', maxX = 120, labelMixSeqRatio = 1 )

mtext('A.',side=3, adj=0, line=1, cex=0.9) #side=1b,2l,3t,4r

#plot margins
par(mar=c(1.5, 0, 2, 0)) #b,l,t,r default c(5, 4, 4, 2)

#effectiveness1 0.8
runcurtis_f2( max_gen=500,  P_1 = 0.01 , P_2 = 0.01 , h.RS1_A0 = 0.5 , h.RS2_0B = 0.5 , exposure = 0.5 , phi.SS1_A0 = 0.8 , phi.SS2_0B = 0.5 , rr_restoration_ins1 = 0.5 , rr_restoration_ins2 = 0.5 , addCombinedStrategy = FALSE, strategyLabels = c('s','','a','m'), ylab="", ylabs = FALSE, xlabs = FALSE, cex.axis = 0.8, addLegend=FALSE, main='', maxX = 120, labelMixSeqRatio = 1 )

mtext('B.',side=3, adj=0, line=1, cex=0.9) #side=1b,2l,3t,4r

#label
par(mar=c(0, 0, 0, 0), pty='m') #b,l,t,r default c(5, 4, 4, 2)
plot( 0, type="n", axes=FALSE, ann=FALSE, pty='m')
mtext('exposure = 0.8', line= -2, side=2, cex=0.9) #side=1b,2l,3t,4r

#plot margins
par(mar=c(1.5, 4, 2, 0)) #b,l,t,r default c(5, 4, 4, 2)

#exposure 0.8
runcurtis_f2( max_gen=500,  P_1 = 0.01 , P_2 = 0.01 , h.RS1_A0 = 0.5 , h.RS2_0B = 0.5 , exposure = 0.8 , phi.SS1_A0 = 0.5 , phi.SS2_0B = 0.5 , rr_restoration_ins1 = 0.5 , rr_restoration_ins2 = 0.5 , addCombinedStrategy = FALSE, strategyLabels = c('s','','a','m'), cex.axis = 0.8, addLegend=FALSE, main='', maxX = 120, labelMixSeqRatio = 1 )

mtext('C.',side=3, adj=0, line=1, cex=0.9) #side=1b,2l,3t,4r


#plot margins
par(mar=c(1.5, 0, 2, 0)) #b,l,t,r default c(5, 4, 4, 2)

#eff1  and  expos 0.8
runcurtis_f2( max_gen=500,  P_1 = 0.01 , P_2 = 0.01 , h.RS1_A0 = 0.5 , h.RS2_0B = 0.5 , exposure = 0.8 , phi.SS1_A0 = 0.8 , phi.SS2_0B = 0.5 , rr_restoration_ins1 = 0.5 , rr_restoration_ins2 = 0.5 , addCombinedStrategy = FALSE, strategyLabels = c('s','','a','m'), ylab="", ylabs = FALSE, cex.axis = 0.8, addLegend=FALSE, main='', maxX = 120, labelMixSeqRatio = 1 )

mtext('D.',side=3, adj=0, line=1, cex=0.9) #side=1b,2l,3t,4r


#adding legend in lower panel
#plot margins
par(mar=c(0, 0, 0, 0), pty='m') #b,l,t,r default c(5, 4, 4, 2)

plot( 0, type="n", axes=FALSE, ann=FALSE, pty='m') 

# legend( 'bottom', legend=c("insecticide1 alone   ", "insecticide1 in mix   ", "insecticide2 alone   ", "insecticide2 in mix   "), col=c("red","red","blue","blue"), lty=c(2,1,2,1), pch=c(NA,NA,NA,NA), bty="n", cex=0.7, horiz =TRUE )

#add for when insecticides same
legend( 'bottom', legend=c("insecticide1 alone   ", "insecticide2 alone   ", "insecticide1 in mix   ", "insecticide2 in mix   ", "insecticides1&2 in mix   "), col=c("red","blue","red","blue","purple"), lty=c(2,2,1,1,1), pch=c(NA,NA,NA,NA,NA), bty="n", cex=0.7, horiz =TRUE )

#label
par(mar=c(0, 0, 0, 0), pty='m') #b,l,t,r default c(5, 4, 4, 2)
plot( 0, type="n", axes=FALSE, ann=FALSE, pty='m')
mtext('effectiveness1 = 0.5', line= -2, side=1, cex=0.9) #side=1b,2l,3t,4r

#label
par(mar=c(0, 0, 0, 0), pty='m') #b,l,t,r default c(5, 4, 4, 2)
plot( 0, type="n", axes=FALSE, ann=FALSE, pty='m')
mtext('effectiveness1 = 0.8', line= -2, side=1, cex=0.9) #side=1b,2l,3t,4r

#4 panel figure looking at effectiveness changes alone


#test of rearranging figures on a grid
#eff1  0.8  1
#eff2
#0.5    1   2
#0.8    3   4


# add side panels for labels
layout( matrix(c(1:6,0,7,7,0,8,9),4,3, byrow = TRUE), heights=c(1,1,0.15,0.2), widths=c(0.1,1,1) )


#label
par(mar=c(0, 0, 0, 0), pty='m') #b,l,t,r default c(5, 4, 4, 2)
plot( 0, type="n", axes=FALSE, ann=FALSE, pty='m')
mtext('effectiveness2 = 0.5', line= -2, side=2, cex=0.9) #side=1b,2l,3t,4r

#base scenario

par(mar=c(1.5, 4, 2, 0)) #b,l,t,r default c(5, 4, 4, 2)

#eff1 0.8
runcurtis_f2( max_gen=500,  P_1 = 0.01 , P_2 = 0.01 , h.RS1_A0 = 0.5 , h.RS2_0B = 0.5 , exposure = 0.5 , phi.SS1_A0 = 0.8 , phi.SS2_0B = 0.5 , rr_restoration_ins1 = 0.5 , rr_restoration_ins2 = 0.5 , addCombinedStrategy = FALSE, strategyLabels = c('s','','a','m'), xlabs = FALSE, cex.axis = 0.8, addLegend=FALSE, main='', maxX = 150, labelMixSeqRatio = 1 )

mtext('A.',side=3, adj=0, line=1, cex=0.9) #side=1b,2l,3t,4r

#plot margins
par(mar=c(1.5, 0, 2, 0)) #b,l,t,r default c(5, 4, 4, 2)

#eff1 1 eff2 0.5
runcurtis_f2( max_gen=500,  P_1 = 0.01 , P_2 = 0.01 , h.RS1_A0 = 0.5 , h.RS2_0B = 0.5 , exposure = 0.5 , phi.SS1_A0 = 1 , phi.SS2_0B = 0.5 , rr_restoration_ins1 = 0.5 , rr_restoration_ins2 = 0.5 , addCombinedStrategy = FALSE, strategyLabels = c('s','','a','m'), ylab="", ylabs = FALSE, xlabs = FALSE, cex.axis = 0.8, addLegend=FALSE, main='', maxX = 150, labelMixSeqRatio = 1 )

mtext('B.',side=3, adj=0, line=1, cex=0.9) #side=1b,2l,3t,4r

#label
par(mar=c(0, 0, 0, 0), pty='m') #b,l,t,r default c(5, 4, 4, 2)
plot( 0, type="n", axes=FALSE, ann=FALSE, pty='m')
mtext('effectiveness2 = 0.8', line= -2, side=2, cex=0.9) #side=1b,2l,3t,4r

#plot margins
par(mar=c(1.5, 4, 2, 0)) #b,l,t,r default c(5, 4, 4, 2)

#eff1 0.8 eff2 0.8
runcurtis_f2( max_gen=500,  P_1 = 0.01 , P_2 = 0.01 , h.RS1_A0 = 0.5 , h.RS2_0B = 0.5 , exposure = 0.5 , phi.SS1_A0 = 0.8 , phi.SS2_0B = 0.8 , rr_restoration_ins1 = 0.5 , rr_restoration_ins2 = 0.5 , addCombinedStrategy = FALSE, strategyLabels = c('s','','a','m'), cex.axis = 0.8, addLegend=FALSE, main='', maxX = 150, labelMixSeqRatio = 1 )

mtext('C.',side=3, adj=0, line=1, cex=0.9) #side=1b,2l,3t,4r


#plot margins
par(mar=c(1.5, 0, 2, 0)) #b,l,t,r default c(5, 4, 4, 2)

#eff1 1 eff2 0.8
runcurtis_f2( max_gen=500,  P_1 = 0.01 , P_2 = 0.01 , h.RS1_A0 = 0.5 , h.RS2_0B = 0.5 , exposure = 0.5 , phi.SS1_A0 = 1 , phi.SS2_0B = 0.8 , rr_restoration_ins1 = 0.5 , rr_restoration_ins2 = 0.5 , addCombinedStrategy = FALSE, strategyLabels = c('s','','a','m'), ylab="", ylabs = FALSE, cex.axis = 0.8, addLegend=FALSE, main='', maxX = 150, labelMixSeqRatio = 1 )

mtext('D.',side=3, adj=0, line=1, cex=0.9) #side=1b,2l,3t,4r


#adding legend in lower panel
#plot margins
par(mar=c(0, 0, 0, 0), pty='m') #b,l,t,r default c(5, 4, 4, 2)

plot( 0, type="n", axes=FALSE, ann=FALSE, pty='m') 

# legend( 'bottom', legend=c("insecticide1 alone   ", "insecticide1 in mix   ", "insecticide2 alone   ", "insecticide2 in mix   "), col=c("red","red","blue","blue"), lty=c(2,1,2,1), pch=c(NA,NA,NA,NA), bty="n", cex=0.7, horiz =TRUE )

#add for when insecticides same
legend( 'bottom', legend=c("insecticide1 alone   ", "insecticide2 alone   ", "insecticide1 in mix   ", "insecticide2 in mix   ", "insecticides1&2 in mix   "), col=c("red","blue","red","blue","purple"), lty=c(2,2,1,1,1), pch=c(NA,NA,NA,NA,NA), bty="n", cex=0.7, horiz =TRUE )

#label
par(mar=c(0, 0, 0, 0), pty='m') #b,l,t,r default c(5, 4, 4, 2)
plot( 0, type="n", axes=FALSE, ann=FALSE, pty='m')
mtext('effectiveness1 = 0.8', line= -2, side=1, cex=0.9) #side=1b,2l,3t,4r

#label
par(mar=c(0, 0, 0, 0), pty='m') #b,l,t,r default c(5, 4, 4, 2)
plot( 0, type="n", axes=FALSE, ann=FALSE, pty='m')
mtext('effectiveness1 = 1', line= -2, side=1, cex=0.9) #side=1b,2l,3t,4r

#4 panel figure looking at dominance and resistance restoration

#rr1  0.5  0.8
#dom1
#0.5    1   2
#0.8    3   4

# add side panels for labels
layout( matrix(c(1:6,0,7,7,0,8,9),4,3, byrow = TRUE), heights=c(1,1,0.15,0.2), widths=c(0.1,1,1) )


#testing
# plot(1)
# plot(1)
# plot(1)
# plot(1)


#label
par(mar=c(0, 0, 0, 0), pty='m') #b,l,t,r default c(5, 4, 4, 2)
plot( 0, type="n", axes=FALSE, ann=FALSE, pty='m')
mtext('dominance1 = 0.5', line= -2, side=2, cex=0.9) #side=1b,2l,3t,4r

#base scenario

par(mar=c(1.5, 4, 2, 0)) #b,l,t,r default c(5, 4, 4, 2)

runcurtis_f2( max_gen=500,  P_1 = 0.01 , P_2 = 0.01 , h.RS1_A0 = 0.5 , h.RS2_0B = 0.5 , exposure = 0.5 , phi.SS1_A0 = 0.5 , phi.SS2_0B = 0.5 , rr_restoration_ins1 = 0.5 , rr_restoration_ins2 = 0.5 , addCombinedStrategy = FALSE, strategyLabels = c('s','','a','m'), xlabs = FALSE, cex.axis = 0.8, addLegend=FALSE, main='', maxX = 120, labelMixSeqRatio = 1 )

mtext('A.',side=3, adj=0, line=1, cex=0.9) #side=1b,2l,3t,4r

#plot margins
par(mar=c(1.5, 0, 2, 0)) #b,l,t,r default c(5, 4, 4, 2)

#rr1 0.8
runcurtis_f2( max_gen=500,  P_1 = 0.01 , P_2 = 0.01 , h.RS1_A0 = 0.5 , h.RS2_0B = 0.5 , exposure = 0.5 , phi.SS1_A0 = 0.5 , phi.SS2_0B = 0.5 , rr_restoration_ins1 = 0.8 , rr_restoration_ins2 = 0.5 , addCombinedStrategy = FALSE, strategyLabels = c('s','','a','m'), ylab="", ylabs = FALSE, xlabs = FALSE, cex.axis = 0.8, addLegend=FALSE, main='', maxX = 120, labelMixSeqRatio = 1 )

mtext('B.',side=3, adj=0, line=1, cex=0.9) #side=1b,2l,3t,4r

#label
par(mar=c(0, 0, 0, 0), pty='m') #b,l,t,r default c(5, 4, 4, 2)
plot( 0, type="n", axes=FALSE, ann=FALSE, pty='m')
mtext('dominance1 = 0.8', line= -2, side=2, cex=0.9) #side=1b,2l,3t,4r

#plot margins
par(mar=c(1.5, 4, 2, 0)) #b,l,t,r default c(5, 4, 4, 2)

#dominance1 0.8
runcurtis_f2( max_gen=500,  P_1 = 0.01 , P_2 = 0.01 , h.RS1_A0 = 0.8 , h.RS2_0B = 0.5 , exposure = 0.5 , phi.SS1_A0 = 0.5 , phi.SS2_0B = 0.5 , rr_restoration_ins1 = 0.5 , rr_restoration_ins2 = 0.5 , addCombinedStrategy = FALSE, strategyLabels = c('s','','a','m'), cex.axis = 0.8, addLegend=FALSE, main='', maxX = 120, labelMixSeqRatio = 1 )

mtext('C.',side=3, adj=0, line=1, cex=0.9) #side=1b,2l,3t,4r


#plot margins
par(mar=c(1.5, 0, 2, 0)) #b,l,t,r default c(5, 4, 4, 2)

#dom1 & rr1 0.8
runcurtis_f2( max_gen=500,  P_1 = 0.01 , P_2 = 0.01 , h.RS1_A0 = 0.8 , h.RS2_0B = 0.5 , exposure = 0.5 , phi.SS1_A0 = 0.5 , phi.SS2_0B = 0.5 , rr_restoration_ins1 = 0.8 , rr_restoration_ins2 = 0.5 , addCombinedStrategy = FALSE, strategyLabels = c('s','','a','m'), ylab="", ylabs = FALSE, cex.axis = 0.8, addLegend=FALSE, main='', maxX = 120, labelMixSeqRatio = 1 )

mtext('D.',side=3, adj=0, line=1, cex=0.9) #side=1b,2l,3t,4r


#adding legend in lower panel
#plot margins
par(mar=c(0, 0, 0, 0), pty='m') #b,l,t,r default c(5, 4, 4, 2)

plot( 0, type="n", axes=FALSE, ann=FALSE, pty='m') 

# legend( 'bottom', legend=c("insecticide1 alone   ", "insecticide1 in mix   ", "insecticide2 alone   ", "insecticide2 in mix   "), col=c("red","red","blue","blue"), lty=c(2,1,2,1), pch=c(NA,NA,NA,NA), bty="n", cex=0.7, horiz =TRUE )

#add for when insecticides same
legend( 'bottom', legend=c("insecticide1 alone   ", "insecticide2 alone   ", "insecticide1 in mix   ", "insecticide2 in mix   ", "insecticides1&2 in mix   "), col=c("red","blue","red","blue","purple"), lty=c(2,2,1,1,1), pch=c(NA,NA,NA,NA,NA), bty="n", cex=0.7, horiz =TRUE )

#label
par(mar=c(0, 0, 0, 0), pty='m') #b,l,t,r default c(5, 4, 4, 2)
plot( 0, type="n", axes=FALSE, ann=FALSE, pty='m')
mtext('resistance restoration1 = 0.5', line= -2, side=1, cex=0.9) #side=1b,2l,3t,4r

#label
par(mar=c(0, 0, 0, 0), pty='m') #b,l,t,r default c(5, 4, 4, 2)
plot( 0, type="n", axes=FALSE, ann=FALSE, pty='m')
mtext('resistance restoration1 = 0.8', line= -2, side=1, cex=0.9) #side=1b,2l,3t,4r

#4 panel figure looking at effectiveness and starting freq.

#eff1  0.5  0.8
#strt freq1
#0.01    1   2
#0.001    3   4

# add side panels for labels
layout( matrix(c(1:6,0,7,7,0,8,9),4,3, byrow = TRUE), heights=c(1,1,0.15,0.2), widths=c(0.1,1,1) )


#label
par(mar=c(0, 0, 0, 0), pty='m') #b,l,t,r default c(5, 4, 4, 2)
plot( 0, type="n", axes=FALSE, ann=FALSE, pty='m')
mtext('frequency1 = 1%', line= -2, side=2, cex=0.9) #side=1b,2l,3t,4r

#base scenario

par(mar=c(1.5, 4, 2, 0)) #b,l,t,r default c(5, 4, 4, 2)

runcurtis_f2( max_gen=500,  P_1 = 0.01 , P_2 = 0.01 , h.RS1_A0 = 0.5 , h.RS2_0B = 0.5 , exposure = 0.5 , phi.SS1_A0 = 0.5 , phi.SS2_0B = 0.5 , rr_restoration_ins1 = 0.5 , rr_restoration_ins2 = 0.5 , addCombinedStrategy = FALSE, strategyLabels = c('s','','a','m'), xlabs = FALSE, cex.axis = 0.8, addLegend=FALSE, main='', maxX = 150, labelMixSeqRatio = 1 )

mtext('A.',side=3, adj=0, line=1, cex=0.9) #side=1b,2l,3t,4r

#plot margins
par(mar=c(1.5, 0, 2, 0)) #b,l,t,r default c(5, 4, 4, 2)

#eff1 0.8
runcurtis_f2( max_gen=500,  P_1 = 0.01 , P_2 = 0.01 , h.RS1_A0 = 0.5 , h.RS2_0B = 0.5 , exposure = 0.5 , phi.SS1_A0 = 0.8 , phi.SS2_0B = 0.5 , rr_restoration_ins1 = 0.5 , rr_restoration_ins2 = 0.5 , addCombinedStrategy = FALSE, strategyLabels = c('s','','a','m'), ylab="", ylabs = FALSE, xlabs = FALSE, cex.axis = 0.8, addLegend=FALSE, main='', maxX = 150, labelMixSeqRatio = 1 )

mtext('B.',side=3, adj=0, line=1, cex=0.9) #side=1b,2l,3t,4r

#label
par(mar=c(0, 0, 0, 0), pty='m') #b,l,t,r default c(5, 4, 4, 2)
plot( 0, type="n", axes=FALSE, ann=FALSE, pty='m')
mtext('frequency1 = 0.1%', line= -2, side=2, cex=0.9) #side=1b,2l,3t,4r

#plot margins
par(mar=c(1.5, 4, 2, 0)) #b,l,t,r default c(5, 4, 4, 2)

#freq1 0.001
runcurtis_f2( max_gen=500,  P_1 = 0.001 , P_2 = 0.01 , h.RS1_A0 = 0.5 , h.RS2_0B = 0.5 , exposure = 0.5 , phi.SS1_A0 = 0.5 , phi.SS2_0B = 0.5 , rr_restoration_ins1 = 0.5 , rr_restoration_ins2 = 0.5 , addCombinedStrategy = FALSE, strategyLabels = c('s','','a','m'), cex.axis = 0.8, addLegend=FALSE, main='', maxX = 150, labelMixSeqRatio = 1 )

mtext('C.',side=3, adj=0, line=1, cex=0.9) #side=1b,2l,3t,4r


#plot margins
par(mar=c(1.5, 0, 2, 0)) #b,l,t,r default c(5, 4, 4, 2)

#freq1 0.001, eff1 0.8
runcurtis_f2( max_gen=500,  P_1 = 0.001 , P_2 = 0.01 , h.RS1_A0 = 0.5 , h.RS2_0B = 0.5 , exposure = 0.5 , phi.SS1_A0 = 0.8 , phi.SS2_0B = 0.5 , rr_restoration_ins1 = 0.5 , rr_restoration_ins2 = 0.5 , addCombinedStrategy = FALSE, strategyLabels = c('s','','a','m'), ylab="", ylabs = FALSE, cex.axis = 0.8, addLegend=FALSE, main='', maxX = 150, labelMixSeqRatio = 1 )

mtext('D.',side=3, adj=0, line=1, cex=0.9) #side=1b,2l,3t,4r


#adding legend in lower panel
#plot margins
par(mar=c(0, 0, 0, 0), pty='m') #b,l,t,r default c(5, 4, 4, 2)

plot( 0, type="n", axes=FALSE, ann=FALSE, pty='m') 

# legend( 'bottom', legend=c("insecticide1 alone   ", "insecticide1 in mix   ", "insecticide2 alone   ", "insecticide2 in mix   "), col=c("red","red","blue","blue"), lty=c(2,1,2,1), pch=c(NA,NA,NA,NA), bty="n", cex=0.7, horiz =TRUE )

#add for when insecticides same
legend( 'bottom', legend=c("insecticide1 alone   ", "insecticide2 alone   ", "insecticide1 in mix   ", "insecticide2 in mix   ", "insecticides1&2 in mix   "), col=c("red","blue","red","blue","purple"), lty=c(2,2,1,1,1), pch=c(NA,NA,NA,NA,NA), bty="n", cex=0.7, horiz =TRUE )

#label
par(mar=c(0, 0, 0, 0), pty='m') #b,l,t,r default c(5, 4, 4, 2)
plot( 0, type="n", axes=FALSE, ann=FALSE, pty='m')
mtext('effectiveness1 = 0.5', line= -2, side=1, cex=0.9) #side=1b,2l,3t,4r

#label
par(mar=c(0, 0, 0, 0), pty='m') #b,l,t,r default c(5, 4, 4, 2)
plot( 0, type="n", axes=FALSE, ann=FALSE, pty='m')
mtext('effectiveness1 = 0.8', line= -2, side=1, cex=0.9) #side=1b,2l,3t,4r

#4 panel figure looking at effectiveness and exposure


#test of rearranging figures on a grid
#eff1  0.5  0.8
#exp
#0.5    1   2
#0.8    3   4

#layout(matrix(1:4,2,2, byrow = TRUE))

#add lower panel for legend
#layout(matrix(c(1,2,3,4,5,5),3,2, byrow = TRUE), heights=c(1,1,0.5), respect=FALSE)

#layout( matrix(c(1,2,3,4,5,5),3,2, byrow = TRUE), heights=c(1,1,0.2), 
# respect=matrix(c(1,1,1,1,0,0),3,2, byrow = TRUE) )
#layout.show(n=5)

# add side panels for labels
layout( matrix(c(1:6,0,7,7,0,8,9),4,3, byrow = TRUE), heights=c(1,1,0.15,0.2), widths=c(0.1,1,1) )

#effectiveness on mixture


#label
par(mar=c(0, 0, 0, 0), pty='m') #b,l,t,r default c(5, 4, 4, 2)
plot( 0, type="n", axes=FALSE, ann=FALSE, pty='m')
mtext('cost 1 = 0', line= -2, side=2, cex=0.9) #side=1b,2l,3t,4r

#base scenario

par(mar=c(1.5, 4, 2, 0)) #b,l,t,r default c(5, 4, 4, 2)

runcurtis_f2( max_gen=500,  P_1 = 0.01 , P_2 = 0.01 , h.RS1_A0 = 0.5 , h.RS2_0B = 0.5 , exposure = 0.5 , phi.SS1_A0 = 0.5 , phi.SS2_0B = 0.5 , rr_restoration_ins1 = 0.5 , rr_restoration_ins2 = 0.5 , addCombinedStrategy = FALSE, strategyLabels = c('s','','a','m'), xlabs = FALSE, cex.axis = 0.8, addLegend=FALSE, main='', maxX = 150, labelMixSeqRatio = 1 )

mtext('A.',side=3, adj=0, line=1, cex=0.9) #side=1b,2l,3t,4r

#plot margins
par(mar=c(1.5, 0, 2, 0)) #b,l,t,r default c(5, 4, 4, 2)

#effectiveness1 0.8
runcurtis_f2( max_gen=500,  P_1 = 0.01 , P_2 = 0.01 , h.RS1_A0 = 0.5 , h.RS2_0B = 0.5 , exposure = 0.5 , phi.SS1_A0 = 0.8 , phi.SS2_0B = 0.5 , rr_restoration_ins1 = 0.5 , rr_restoration_ins2 = 0.5 , addCombinedStrategy = FALSE, strategyLabels = c('s','','a','m'), ylab="", ylabs = FALSE, xlabs = FALSE, cex.axis = 0.8, addLegend=FALSE, main='', maxX = 150, labelMixSeqRatio = 1 )

mtext('B.',side=3, adj=0, line=1, cex=0.9) #side=1b,2l,3t,4r

#label
par(mar=c(0, 0, 0, 0), pty='m') #b,l,t,r default c(5, 4, 4, 2)
plot( 0, type="n", axes=FALSE, ann=FALSE, pty='m')
mtext('cost 1 = 0.15', line= -2, side=2, cex=0.9) #side=1b,2l,3t,4r

#plot margins
par(mar=c(1.5, 4, 2, 0)) #b,l,t,r default c(5, 4, 4, 2)

#exposure 0.8
runcurtis_f2( max_gen=500,  P_1 = 0.01 , P_2 = 0.01 , h.RS1_A0 = 0.5 , h.RS2_0B = 0.5 , exposure = 0.5 , phi.SS1_A0 = 0.5 , phi.SS2_0B = 0.5 , rr_restoration_ins1 = 0.5 , rr_restoration_ins2 = 0.5 , addCombinedStrategy = FALSE, strategyLabels = c('s','','a','m'), xlabs = TRUE, cex.axis = 0.8, addLegend=FALSE, main='', maxX = 150, labelMixSeqRatio = 1, z.RR1_00 = 0.15 )

mtext('C.',side=3, adj=0, line=1, cex=0.9) #side=1b,2l,3t,4r


#plot margins
par(mar=c(1.5, 0, 2, 0)) #b,l,t,r default c(5, 4, 4, 2)

#eff1  and  expos 0.8
runcurtis_f2( max_gen=500,  P_1 = 0.01 , P_2 = 0.01 , h.RS1_A0 = 0.5 , h.RS2_0B = 0.5 , exposure = 0.5 , phi.SS1_A0 = 0.8 , phi.SS2_0B = 0.5 , rr_restoration_ins1 = 0.5 , rr_restoration_ins2 = 0.5 , addCombinedStrategy = FALSE, strategyLabels = c('s','','a','m'), ylab="", ylabs = FALSE, xlabs = TRUE, cex.axis = 0.8, addLegend=FALSE, main='', maxX = 150, labelMixSeqRatio = 1, z.RR1_00 = 0.15  )

mtext('D.',side=3, adj=0, line=1, cex=0.9) #side=1b,2l,3t,4r


#adding legend in lower panel
#plot margins
par(mar=c(0, 0, 0, 0), pty='m') #b,l,t,r default c(5, 4, 4, 2)

plot( 0, type="n", axes=FALSE, ann=FALSE, pty='m') 

# legend( 'bottom', legend=c("insecticide1 alone   ", "insecticide1 in mix   ", "insecticide2 alone   ", "insecticide2 in mix   "), col=c("red","red","blue","blue"), lty=c(2,1,2,1), pch=c(NA,NA,NA,NA), bty="n", cex=0.7, horiz =TRUE )

#add for when insecticides same
legend( 'bottom', legend=c("insecticide1 alone   ", "insecticide2 alone   ", "insecticide1 in mix   ", "insecticide2 in mix   ", "insecticides1&2 in mix   "), col=c("red","blue","red","blue","purple"), lty=c(2,2,1,1,1), pch=c(NA,NA,NA,NA,NA), bty="n", cex=0.7, horiz =TRUE )

#label
par(mar=c(0, 0, 0, 0), pty='m') #b,l,t,r default c(5, 4, 4, 2)
plot( 0, type="n", axes=FALSE, ann=FALSE, pty='m')
mtext('effectiveness1 = 0.5', line= -2, side=1, cex=0.9) #side=1b,2l,3t,4r

#label
par(mar=c(0, 0, 0, 0), pty='m') #b,l,t,r default c(5, 4, 4, 2)
plot( 0, type="n", axes=FALSE, ann=FALSE, pty='m')
mtext('effectiveness1 = 0.8', line= -2, side=1, cex=0.9) #side=1b,2l,3t,4r

# add side panels for labels
#layout( matrix(c(1:6,0,7,7,0,8,9),4,3, byrow = TRUE), heights=c(1,1,0.15,0.2), widths=c(0.1,1,1) )
layout( matrix(c(1:2), byrow = TRUE), heights=c(1,0.15) )

#plot margins
par(mar=c(1.5, 0, 2, 0)) #b,l,t,r default c(5, 4, 4, 2)

# runcurtis_f2( max_gen=500,  P_1 = 0.01 , P_2 = 0.01 , h.RS1_A0 = 0.3 , h.RS2_0B = 0.66 , exposure = 0.5 , phi.SS1_A0 = 0.43 , phi.SS2_0B = 0.98 , rr_restoration_ins1 = 0.63 , rr_restoration_ins2 = 0.84 , addCombinedStrategy = FALSE, strategyLabels = c('s','','a','m'), cex.axis = 0.8, addLegend=FALSE, labelMixSeqRatio = 1 ) 

# Alpha-cypermethrin (pyr) incl cost, exposure 0.8
runcurtis_f2( max_gen=500,  P_1 = 0.01 , P_2 = 0.01, z.RR1_00 = 0.19, h.RS1_00 = 1,  h.RS1_A0 = 0.54 , h.RS2_0B = 0.66 , exposure = 0.8 , phi.SS1_A0 = 0.29 , phi.SS2_0B = 0.98 , rr_restoration_ins1 = 1 , rr_restoration_ins2 = 0.84 , addCombinedStrategy = FALSE, strategyLabels = c('s','','a','m'), cex.axis = 0.8, addLegend=FALSE, labelMixSeqRatio = 1 ) 

#adding legend in lower panel
#plot margins
par(mar=c(0, 0, 0, 0), pty='m') #b,l,t,r default c(5, 4, 4, 2)

plot( 0, type="n", axes=FALSE, ann=FALSE, pty='m') 

legend( 'bottom', legend=c("pyrethroid alone   ", "pyrethroid in mix   ", "carbamate alone   ", "carbamate in mix   "), col=c("red","red","blue","blue"), lty=c(2,1,2,1), pch=c(NA,NA,NA,NA), bty="n", cex=0.7, horiz =TRUE )

# not evaluated at the moment
#4 panel figure looking at male exposure prop and starting freq.

#eff1  0.5  0.8
#maleprop
#0.8    1   2
#0.2    3   4

# add side panels for labels
layout( matrix(c(1:6,0,7,7,0,8,9),4,3, byrow = TRUE), heights=c(1,1,0.15,0.2), widths=c(0.1,1,1) )


#label
par(mar=c(0, 0, 0, 0), pty='m') #b,l,t,r default c(5, 4, 4, 2)
plot( 0, type="n", axes=FALSE, ann=FALSE, pty='m')
mtext('male exposure = 0.8', line= -2, side=2, cex=0.9) #side=1b,2l,3t,4r

#base scenario

par(mar=c(1.5, 4, 2, 0)) #b,l,t,r default c(5, 4, 4, 2)

runcurtis_f2( male_exposure_prop=0.8, max_gen=500,  P_1 = 0.01 , P_2 = 0.01 , h.RS1_A0 = 0.5 , h.RS2_0B = 0.5 , exposure = 0.5 , phi.SS1_A0 = 0.5 , phi.SS2_0B = 0.5 , rr_restoration_ins1 = 0.5 , rr_restoration_ins2 = 0.5 , addCombinedStrategy = FALSE, strategyLabels = c('s','','a','m'), xlabs = FALSE, cex.axis = 0.8, addLegend=FALSE, main='', maxX = 200, labelMixSeqRatio = 1 )

mtext('A.',side=3, adj=0, line=1, cex=0.9) #side=1b,2l,3t,4r

#plot margins
par(mar=c(1.5, 0, 2, 0)) #b,l,t,r default c(5, 4, 4, 2)

#eff1 0.8
runcurtis_f2( male_exposure_prop=0.8, max_gen=500,  P_1 = 0.01 , P_2 = 0.01 , h.RS1_A0 = 0.5 , h.RS2_0B = 0.5 , exposure = 0.5 , phi.SS1_A0 = 0.8 , phi.SS2_0B = 0.5 , rr_restoration_ins1 = 0.5 , rr_restoration_ins2 = 0.5 , addCombinedStrategy = FALSE, strategyLabels = c('s','','a','m'), ylab="", ylabs = FALSE, xlabs = FALSE, cex.axis = 0.8, addLegend=FALSE, main='', maxX = 200, labelMixSeqRatio = 1 )

mtext('B.',side=3, adj=0, line=1, cex=0.9) #side=1b,2l,3t,4r

#label
par(mar=c(0, 0, 0, 0), pty='m') #b,l,t,r default c(5, 4, 4, 2)
plot( 0, type="n", axes=FALSE, ann=FALSE, pty='m')
mtext('male exposure =0.2', line= -2, side=2, cex=0.9) #side=1b,2l,3t,4r

#plot margins
par(mar=c(1.5, 4, 2, 0)) #b,l,t,r default c(5, 4, 4, 2)

#freq1 0.001
runcurtis_f2( male_exposure_prop=0.2, max_gen=500,  P_1 = 0.01 , P_2 = 0.01 , h.RS1_A0 = 0.5 , h.RS2_0B = 0.5 , exposure = 0.5 , phi.SS1_A0 = 0.5 , phi.SS2_0B = 0.5 , rr_restoration_ins1 = 0.5 , rr_restoration_ins2 = 0.5 , addCombinedStrategy = FALSE, strategyLabels = c('s','','a','m'), cex.axis = 0.8, addLegend=FALSE, main='', maxX = 200, labelMixSeqRatio = 1 )

mtext('C.',side=3, adj=0, line=1, cex=0.9) #side=1b,2l,3t,4r


#plot margins
par(mar=c(1.5, 0, 2, 0)) #b,l,t,r default c(5, 4, 4, 2)

#freq1 0.001, eff1 0.8
runcurtis_f2( male_exposure_prop=0.2, max_gen=500,  P_1 = 0.01 , P_2 = 0.01 , h.RS1_A0 = 0.5 , h.RS2_0B = 0.5 , exposure = 0.5 , phi.SS1_A0 = 0.8 , phi.SS2_0B = 0.5 , rr_restoration_ins1 = 0.5 , rr_restoration_ins2 = 0.5 , addCombinedStrategy = FALSE, strategyLabels = c('s','','a','m'), ylab="", ylabs = FALSE, cex.axis = 0.8, addLegend=FALSE, main='', maxX = 200, labelMixSeqRatio = 1 )

mtext('D.',side=3, adj=0, line=1, cex=0.9) #side=1b,2l,3t,4r


#adding legend in lower panel
#plot margins
par(mar=c(0, 0, 0, 0), pty='m') #b,l,t,r default c(5, 4, 4, 2)

plot( 0, type="n", axes=FALSE, ann=FALSE, pty='m') 

# legend( 'bottom', legend=c("insecticide1 alone   ", "insecticide1 in mix   ", "insecticide2 alone   ", "insecticide2 in mix   "), col=c("red","red","blue","blue"), lty=c(2,1,2,1), pch=c(NA,NA,NA,NA), bty="n", cex=0.7, horiz =TRUE )

#add for when insecticides same
legend( 'bottom', legend=c("insecticide1 alone   ", "insecticide2 alone   ", "insecticide1 in mix   ", "insecticide2 in mix   ", "insecticides1&2 in mix   "), col=c("red","blue","red","blue","purple"), lty=c(2,2,1,1,1), pch=c(NA,NA,NA,NA,NA), bty="n", cex=0.7, horiz =TRUE )

#label
par(mar=c(0, 0, 0, 0), pty='m') #b,l,t,r default c(5, 4, 4, 2)
plot( 0, type="n", axes=FALSE, ann=FALSE, pty='m')
mtext('effectiveness1 = 0.5', line= -2, side=1, cex=0.9) #side=1b,2l,3t,4r

#label
par(mar=c(0, 0, 0, 0), pty='m') #b,l,t,r default c(5, 4, 4, 2)
plot( 0, type="n", axes=FALSE, ann=FALSE, pty='m')
mtext('effectiveness1 = 0.8', line= -2, side=1, cex=0.9) #side=1b,2l,3t,4r

\pagebreak

References

AndySouth/resistance documentation built on Nov. 12, 2020, 3:39 a.m.

rdrr.io home R language documentation Run R code online

CRAN packages Bioconductor packages R-Forge packages GitHub packages

Note that we can't provide technical support on individual packages. You should contact the package authors for that.

AndySouth/resistance
an insecticide resistance population genetics model with 2 loci and 2 insecticides

In AndySouth/resistance: an insecticide resistance population genetics model with 2 loci and 2 insecticides

Abstract

Background

Results

Conclusions

Keywords

Background

Methods

Results

Single insecticide

Two insecticides

Insecticide exposure and effectiveness

Resistance restoration and dominance

Starting frequencies of resistance

Costs of resistance

Linkage disequilibrium

Setting model inputs from field data

Discussion

Two insecticides

Difference between exposure and effectiveness

Costs of resistance

Using field derived inputs

Consistency with previous work

Caveats

Poor implementation

Conclusions

Declarations

Ethics approval and consent to participate

Consent for publication

Availability of data and material

Competing interests

Funding

Authors' contributions

Acknowledgements

Figures

References

R Package Documentation

Browse R Packages

We want your feedback!

AndySouth/resistance an insecticide resistance population genetics model with 2 loci and 2 insecticides

In AndySouth/resistance: an insecticide resistance population genetics model with 2 loci and 2 insecticides

Abstract

Background

Results

Conclusions

Keywords

Background

Methods

Results

Single insecticide

Two insecticides

Insecticide exposure and effectiveness

Resistance restoration and dominance

Starting frequencies of resistance

Costs of resistance

Linkage disequilibrium

Setting model inputs from field data

Discussion

Two insecticides

Difference between exposure and effectiveness

Costs of resistance

Using field derived inputs

Consistency with previous work

Caveats

Poor implementation

Conclusions

Declarations

Ethics approval and consent to participate

Consent for publication

Availability of data and material

Competing interests

Funding

Authors' contributions

Acknowledgements

Figures

References

R Package Documentation

Browse R Packages

We want your feedback!

AndySouth/resistance
an insecticide resistance population genetics model with 2 loci and 2 insecticides