Insecticide resistance management strategies for public health control of mosquitoes exhibiting polygenic resistance: A comparison of sequences, rotations, and mixtures

Abstract Malaria control uses insecticides to kill Anopheles mosquitoes. Recent successes in malaria control are threatened by increasing levels of insecticide resistance (IR), requiring insecticide resistance management (IRM) strategies to mitigate this problem. Field trials of IRM strategies are usually prohibitively expensive with long timeframes, and mathematical modeling is often used to evaluate alternative options. Previous IRM models in the context of malaria control assumed IR to have a simple (monogenic) basis, whereas in natural populations, IR will often be a complex polygenic trait determined by multiple genetic variants. A quantitative genetics model was developed to model IR as a polygenic trait. The model allows insecticides to be deployed as sequences (continuous deployment until a defined withdrawal threshold, termed “insecticide lifespan”, as indicated by resistance diagnosis in bioassays), rotations (periodic switching of insecticides), or full‐dose mixtures (two insecticides in one formulation). These IRM strategies were compared based on their “strategy lifespan” (capped at 500 generations). Partial rank correlation and generalized linear modeling was used to identify and quantify parameters driving the evolution of resistance. Random forest models were used to identify parameters offering predictive value for decision‐making. Deploying single insecticides as sequences or rotations usually made little overall difference to their “strategy lifespan”, though rotations displayed lower mean and peak resistances. Deploying two insecticides in a full‐dose mixture formulation was found to extend the “strategy lifespan” when compared to deploying each in sequence or rotation. This pattern was observed regardless of the level of cross resistance between the insecticides or the starting level of resistance. Statistical analysis highlighted the importance of insecticide coverage, cross resistance, heritability, and fitness costs for selecting an appropriate IRM strategy. Full‐dose mixtures appear the most promising of the strategies evaluated, with the longest “strategy lifespans”. These conclusions broadly corroborate previous results from monogenic models.


| INTRODUC TI ON
Long-lasting insecticide-treated nets (LLINs) and indoor residual spraying (IRS) play a dominant role in reducing the burden of malaria (Bhatt et al., 2015). The evolution of insecticide resistance (IR) poses a major threat to sustained reductions in malaria transmission and prevalence (e.g., Hemingway et al., 2016). Pyrethroids are the primary insecticide class used on LLINs, but pyrethroid resistance is geographically widespread (Hancock et al., 2020) and is increasing in frequency and intensity .
How entomological outcomes (i.e., greater mosquito survival) translate into impacts on malaria epidemiology is uncertain and complex (Van Hul et al., 2021). However, tools designed to kill mosquitoes or prevent them from blood-feeding, which cease to be capable of doing so because of resistance (Asidi et al., 2012;Irish et al., 2008), are clearly suboptimal.
The Global Plan for Insecticide Resistance Management (GPIRM) was developed to help mitigate the potential impact IR (WHO, 2012).
Compared to agriculture, public health has relatively few insecticide classes available . Insecticides for malaria control are mainly deployed to target adult mosquitoes indoors as either LLINs or IRS, requiring the insecticides to be safe for close human contact, which constrains the development of novel products. LLINs are designed to have 3-year lifespans, whilst IRS is usually deployed annually (WHO, 2020). Therefore, insecticides used in an insecticide resistance management (IRM) strategy for public health cannot be changed rapidly, unlike agriculture where re-application can occur monthly or even weekly. This necessity for safe and longlasting insecticides has led to an over-reliance on a small number of insecticides, predominantly pyrethroids (Oxborough, 2016).
There is a need to evaluate IRM strategies designed to prolong the operational lifespan of insecticides to ensure future insecticides are not rapidly lost to resistance. This is especially important considering the new insecticides being developed by the Innovative Vector Control Consortium and industry partners. Evaluating IRM strategies in the laboratory and field is challenging due to the number of potential IRM strategies, the need for replication over different ecological/epidemiological settings, and the long study durations needed to observe phenotypic changes in resistance. Cluster randomized control trials used for evaluating the epidemiological effectiveness of malaria vector control tools generally last no longer than 2-3 years (Wilson et al., 2015), which is unlikely to be of sufficient duration to detect the benefit of any one IRM strategy over another.
Mathematical modeling and computer simulations can simulate the IR response to IRM strategies over decades, and therefore provide a valuable evaluation tool (Tabashnik, 1986). Computer simulations can be run over numerous scenarios to provide insights into which IRM strategies are optimal/sub-optimal under different ecological, epidemiological, and operational contexts. Evaluating IRM strategies using computer simulations allows for the identification of potential problems with IRM strategies prior to implementation in the field where expensive and time-consuming errors can be made.
Theoretical and mathematical modeling of IRM for agriculture and public health has a long history (e.g., Comins, 1977;Curtis, 1985;Mani, 1985;Wood, 1981). However much of this previous work has modelled IR as a monogenic trait i.e., a trait encoded by a single mutation in a single gene. Few modeling studies have considered IR to be a polygenic trait (but see, for example, Gardner et al., 1998;Haridas & Tenhumberg, 2018;Via, 1986), of which none have evaluated IRM strategies in a public health context.
In more recent years, mathematical models have been developed which take advantage of increased computational power allowing for improved evaluation of IRM strategies in agriculture (e.g., Helps et al., 2017) and public health (e.g., Hastings et al., 2022;Levick et al., 2017;Madgwick & Kanitz, 2022b), of which these examples are monogenic models.
Public health IRM recommendations from models are therefore currently limited to scenarios where resistance in the field is monogenic. However mosquitoes are known to harbor a diverse array of resistance genes and mechanisms, including target site, metabolic, cuticular and behavioral resistance mechanisms (Balabanidou et al., 2018), alongside nonspecific, largely environmentally influenced factors associated with increased tolerance to insecticides, for example mosquito age, and physiological condition (Lissenden et al., 2021).
Resistance management strategies often have confusing naming systems resulting from the often-independent development and evaluation of resistance management strategies across disciplines (Peck, 2001;Rex Consortium, 2007). Therefore, we define our terminology here for the three IRM strategies to investigate, and how this relates to terms used in previous work: • Sequences, involve continuously deploying a particular insecticide formulation until a designated level of resistance is reached (the withdrawal threshold), at which point it is replaced by another insecticide to which there is less resistance. This strategy is also referred to as responsive alternation or series application in the literature (Rex Consortium, 2013). The threshold for withdrawing the current insecticide and replacing with the next insecticide the sequence is usually based on a bioassay survival threshold or, in the case of monogenic traits, an allele frequency threshold.
• Rotations, also referred to as cycling, alternation, or periodic application (Rex Consortium, 2013), involve the periodic preplanned switching between insecticide formulations over time, such that one insecticide is temporarily replaced with another insecticide.
• Mixtures are confusingly also occasionally referred to as combinations (Rex Consortium, 2013) or pyramiding when referring to transgenic crops (Roush, 1998). In public health, the term combination would more generally refer to the use of two different insecticides within the same household, but not in the same control tool, such as deploying both LLINs and IRS (WHO, 2012). We define mixtures as the use of two insecticides in a single formulation such that the target insect inevitably encounters both insecticides simultaneously.
A comprehensive review of IRM modeling highlighted a lack of models including quantitative resistance and cross resistance (Rex Consortium, 2010). Quantitative genetic modeling assumes many genes determine IR level and each gene has only a minor phenotypic effect. It is the accumulated effect of these genes which combine to give larger phenotypic effects (Walsh & Lynch, 2018). In this study we model IR as a quantitative trait, assuming a polygenic basis of resistance and allowing for cross-resistance between insecticides.
We present a flexible quantitative genetics model ("polyres") for the evaluation of IRM strategies. The model is calibrated primarily for Anopheles gambiae but can be adapted for other vector species.
We investigate how insecticides could be deployed temporally, as either sequences, rotations, or mixtures. Based on the results from monogenic models, we would hypothesize the difference between the rotation and sequence strategies to be small, and for the mixture strategy to perform best.

| An overview of the quantitative genetics model
In the "polyres" model, IR is assumed to be a classically quantitative genetic trait encoded by a large number of genes, each with very small effect on the phenotype (this assumption means genetic variances do not change substantially over selection). The model assumes discrete nonoverlapping generations of mosquitoes, a standard assumption in both quantitative and population genetics modeling.
Female mosquitoes lay eggs, which hatch, develop, and mature into adults. Adult mosquitoes of both sexes may then encounter the deployed insecticide(s). Female mosquitoes then lay the eggs that constitute the next discrete generation and then die, that is, the model allows for only a single insecticide exposure per generation.
We assume mating between male and female mosquitoes occurs after the insecticide encounters, that is, after insecticide selection which is consistent with equivalent monogenic models (e.g., Levick et al., 2017;Madgwick & Kanitz, 2022b), and all females successfully mate. This ordering of insecticide exposure, mating, then dispersal is important as this impacts the level of selection (Sudo et al., 2018

| Developing a quantitative IR scale
To develop a polygenic model which allows for the tracking of quantitative IR, the first step is to develop an underlying scale of IR that can be converted to a measurable phenotype, in this case mosquito survival to insecticide exposure (see Supplement S1 for details). In the field, the IR status of a mosquito population is typically measured using standardized diagnostic dose bioassays such as CDC bottle or WHO cylinder bioassays. These bioassays give the proportion of mosquitoes surviving contact with the insecticide (the measurable phenotype). We created a scale termed the Polygenic Resistance Score (PRS), denoted by z, which we use as a measure of resistance, and is described in more detail below. This underlying PRS scale produces a bioassay survival phenotype described by the Hill-variant of the Michaelis-Menten equation (where n = 1) which converts the PRS to bioassay survival.
i is the proportion of mosquitoes surviving in a diagnostic dose bioassay to insecticide i when they have a population mean PRS of z I for corresponding resistance trait I. K max is the maximum proportion of mosquitoes that could survive in a bioassay, which is 1. The z 50 value is the PRS that gives 50% bioassay survival. For ease of interpretability, we have scaled the PRS such that z I =100 equates to 10% bioassay survival. The z 50 for this scale is therefore calculated from Equation 1a to be 900. Our criterion for withdrawing an insecticide is >10% bioassay survival (<90% mortality) which is defined as the point when resistance is confirmed (WHO, 2018) so the evolution of resistance occurs over a scale of z = 0 (the starting point) to z = 100 (the withdrawal threshold).
We note our withdrawal threshold is set optimistically low, and in part this is done to help "force" one strategy to be better than other, if the withdrawal threshold is set high then all strategies would run to the end of the simulation and therefore all appear equally good.
This relationship between PRS and bioassay survival is shown in Figure 1. This calibration of K max = 1 and z 50 = 900 is stable (i.e., bioassay mortality = 10% when z = 100) for population standard deviations up to 25 (Supplement S3, Table S2). A full description of all model parameters can be found in Table 3.

| Converting bioassay survival to field survival
The PRS must then be converted from bioassay survival to field survival. As mosquito survival in bioassays has been found to be correlated with mortality in experimental huts (Churcher et al., 2016), we can update Equation 1a: where K F i is the survival to insecticide i in experimental huts, which we use as an approximation for field survival. K B i is the bioassay survival to There is only dispersal between the intervention site and refugia Standard assumption in two-patch models (see Comins, 1977) 10

| The response to selection
The changes in population mean PRS (z) are tracked using the Breeder's equation. The Breeder's equation is synonymously known as the Lush equation (Walsh & Lynch, 2018): R is the response to selection, that is, the inherited change in the mean value of the PRS between generations. h 2 is the narrow sense heritability of the trait and S is the selection differential. The selection differential is the within-generation change of the mean PRS of the population due to insecticide selection.
where z P is the mean PRS of the parents (i.e., those that either do not encounter insecticide, or encounter and survive insecticide, and go on to produce offspring), and z is the mean PRS prior to selection. Only female mosquitoes blood-feed, and because (2a) R = Sh 2 (2ai) S = z P − z TA B L E 2 Headline summary of polyres model methods.

Methods section
Headline summary

Mean Polygenic Resistance Score in the Refugia
The mean Polygenic Resistance Score of the mosquito population in the refugia before selection has occurred that generation Internally calculated The following input parameters were allowed to vary in the simulations:

ΓI
Degree of Cross-Resistance Degree of cross-resistance from one insecticide on the trait not associated with that insecticide −0.5 to 0.5 at 0.1 intervals h 2 Heritability The heritability of a polygenic trait Uniform 0.05 to 0.30 x Female Insecticide Exposure Proportion of female mosquitoes in the intervention site that encounter and are exposed to the deployed insecticide The selection differential for females and males is: where x is exposure to the insecticide (proportion of females encountering the insecticide in a generation) and beta is a scaling factor (see below). The term m is the proportion of male An. gambiae exposed to the insecticide as a proportion of the female mosquitoes, which therefore accounts for female An. gambiae mosquitoes being more likely to encounter insecticides especially in the form of an LLIN or IRS. The male and female selection pressures can then be implemented back into Equation 2b: Note that insecticide deployments that do not differentially impact males and females (e.g., spraying larval breeding sites) could be investigated simply by setting m = 1.
An important operational difficulty in applying the breeder's equation (i.e., 2a) to selection in the field is that neither the selection differentials imposed by insecticides in the field, nor the field heritability of traits are known (heritability in the field will be much lower than in the laboratory due to the much larger environmental variance in field heritability). We therefore included an exposure scaling factor ( ) to incorporate this uncertainty. Details of the calibration (Figure 2) of are found in Supplement S1.

| Fitness costs
Insecticide resistant mosquitoes may be less fit than insecticide sus- It is assumed that the fitness costs remain constant regardless of the magnitude of the IR level, z. Note that we assume that z = 0 is evolutionary stable (i.e., resistance levels are at their basal, predeployment level) and fitness costs cannot reduce z to less than zero, for example when the insecticide is not being deployed. Additional details are available in Supplement S1.

| Accounting for mosquito dispersal
Mosquitoes may disperse between the intervention site and an insecticide-free refugia. Dispersal from intervention to refugia is given in Equation 6a and from refugia to intervention site Equation 6b (see supplementary Information in Hastings et al., 2022 for details).
where r t is the number of mosquitoes migrating from the insecticidetreated location to the insecticide-free refugia. and r u is the number of mosquitoes migrating from the insecticide-free refugia to the insecticide-treated location. C is the coverage of the insecticide (i.e., the proportion of the population covered by the intervention), and e is the proportion of mosquitoes dispersing.

| Tracking the polygenic resistance score
The impact of the insecticide selection pressure and fitness costs associated with IR during insecticide deployment in the location where the insecticide is deployed is quantified as the change over a generation and is given by the response to insecticide selection and the fitness costs of IR as explained previously.
where the superscript I represents the insecticide-resistance trait (resistance to insecticide i), and the subscript t represents the intervention site.
The proportion of the total mosquito population that is covered by the intervention site When an insecticide is not deployed in the treatment location (e.g., because it has been rotated out) there is no insecticide selection but there are fitness costs which are calculated as explained previously.
The impact of mosquito dispersal is then allowed as in Equation 7b i.e.,: In the insecticide-free refugia (subscript u), only IR costs are present as there is never any insecticide deployment.
Migration is also allowed to occur, with mosquitoes dispersing from or to the refugia.
where z I u ′ 1 − r u is the mean resistance of the individuals staying in the refugia site and z I u ′ r t is the mean resistance of those joining from the intervention site.
Note the use of primes: the single prime (e.g., z I u ′ ) indicates the value of z after insecticide selection and the double prime (e.g., z I u ′′ is the value after both insecticide selection and mosquito dispersal.
It is the double prime value that becomes the mean PRS value in the next generation.

| Special case 1: Cross resistance and crossselection between insecticides
Cross resistance between insecticides is often considered an important factor when selecting insecticides for an IRM strategy especially relevant to the mixture strategy (e.g., Curtis, 1985), yet is often not included in mathematical models evaluating IRM strategies (Rex Consortium, 2010). Via (1986) included genetic correlation, but in a model which tracked the median lethal dose (LD 50 ) to the respective insecticides and therefore needed a more complex method to calculate the genetic correlation. We note the terms cross selection and cross resistance are often used interchangeably, but to aid our readers we use the term cross resistance as we expect readers to be more familiar with this concept.
The degree of genetic correlation between the level of IR to insecticide i (trait I) and (trait Γ) is quantified as IΓ . As all the traits measured in this model are on the same scale (e.g., z = 100 for trait I is 10% bioassay survival to insecticide i , and z = 100 for trait J is 10% bioassay survival to insecticide j), there is no need for regression/ variance coefficients to translate between different scales.
Histogram of the Insecticide Lifespan Assuming the Exposure Scaling Factor (β) = 10. Calibration of the exposure scaling factor is needed to convert selection to the desired timescale of expected insecticide lifespans (i.e., Equation 4). The exposure scaling factor was varied until the peak of the histogram was approximately 10 years, as we expect insecticides to last this length. The red vertical lines are at 8 and 12 years, indicating the range inside which we would expect insecticides to most frequently reached the reach the 10% bioassay survival withdrawal threshold and therefore reach the end of its insecticide lifespan. Equation 8a is updated to include cross resistance/selection, such that when insecticide i is not deployed and insecticide is deployed, there is also selection on trait I from insecticide based on the genetic correlation trait Γ and trait I, termed ΓI .
where ∈ {j, k, l … }, the set of the other insecticides that may be deployed in the mixture.
Where Γ ∈ {J, K, L … }, the set of traits that correspond to IR to the corresponding insecticide.
Equation 8a(i) is passed to Equation 9a as previously described to allow migration to occur as previously described. The refugia equations do not need updating as we assume no insecticides under consideration in the IRM are deployed in the refugia.

| Special case 2: Inclusion of mixtures
Insecticide mixtures describe the simultaneous deployment of two (or, in theory, more) insecticides in the same formulation. Given that LLINs and IRS formulations currently submitted for WHO prequalification do not include more than two insecticides, mixtures are limited to contain only two insecticides in our simulations. The primary idea behind mixtures is that if one insecticide in the mixture fails to kill the mosquito, the other insecticide of the mixture will do so. This requires Equation 7a to be updated. When an insecticide mixture contains both insecticides i and j the mosquito must survive the encounter with one part of the insecticide before selection can occur. In this example, prior to selection on trait I the mosquito must first survive its encounter with insecticide j.
where K F j is the average probability of a mosquito surviving exposure to insecticide j in the field when deployed at the full dose, which depends on the concurrent value of z J and is obtained from Equation 1b. When insecticide i is not deployed in a mixture, then the model proceeds through Equation 8a as normal. The refugia Equations 9a and 9b remained unchanged apart from being passed the updated z I t ′ from Equation 7a(ii).

| Special case 3: Cross resistance and selection with mixtures
When an insecticide mixture is deployed that consists of insecticide i and insecticide j and there is cross resistance between insecticides there is both direct and indirect selection on trait I which is then scaled by the survival probability to the second insecticide in the mixture.

Equation 7a
(ii) is therefore updated so that the survival probability to the second insecticide in the mixture acts on the changes in resistance intensity.
where the second term incorporates direct selection, and the third term incorporates indirect selection by cross resistance. Note that if cross-selection is absent ΓI = 0 so Equation 7a(iii) is, as expected, equivalent to Equation 7a(ii). These are then passed to Equations 7b and 8b to allow for migration. The refugia Equations 9a and 9b, remain unchanged besides being passed the updated z I t ′ from either Equations 7a(iii) or 8a(ii).
The mathematical model described above is coded in R (R Core Team, 2020), version 4.0.3. Model code is written following modern coding practices including the use of modular coding, unit testing and maintenance in a version-controlled repository.

| Description of the simulations
We define several key terms and model rules here: Insecticide Armory: The number of different insecticide formulations available for deployment. Only insecticide formulations in the armory can be deployed. Insecticides can be withdrawn from the armory and returned to the armory based on the withdrawal threshold and return threshold described below.
Withdrawal Threshold: Is the resistance threshold whereupon an insecticide is considered to have failed and is withdrawn from the armory. We set the withdrawal threshold at 10% bioassay survival (90% bioassay mortality) as a ≥10% bioassay survival indicates there is confirmed resistance in the mosquito population (WHO, 2018).
The withdrawal threshold can be set by the user to be higher or lower.
Return Threshold: The bioassay survival an insecticide must reach for a previously failed insecticide to be return to the insecticide armory and become re-available for deployment. This is to prevent an insecticide that has recently failed being immediately returned to deployment. In our presented simulations the return threshold was set at 8% bioassay survival, though this can be set by the user. For example, if an insecticide reaches the withdrawal threshold of 10% bioassay survival the insecticide remains unavailable for re-deployment until its bioassay survival falls below 8%. We use 8% as this value allows withdrawn insecticides to have the potential redeployed, as lower return thresholds are unlikely to ever be reached in the simulations. The withdrawal threshold value can therefore impact the difference between sequences and rotations, if set too low a withdrawn insecticide (in sequence) is not returned to the armory before the second insecticide is withdrawn. An illustrative example of the process of withdrawing and returning insecticides is presented in the sequence simulation example in Figure 3.  Procurement/logistical failures leading to the deployment of just a solo insecticide, and therefore the insecticide is still available.
For simulations where insecticides are given unique properties, we must adapt the rotations strategy so that this strategy does not appear artificially ineffective due to the previously described rotations rules: • Adaptive rotations: In adaptative rotations we relax the restriction that prevents an insecticide being immediately redeployed. This means when only one insecticide is available for deployment (because the other is above the return threshold), it continues to be deployed. However, at the first deployment opportunity where the other insecticide becomes available (i.e., has fallen below the return threshold), the strategy defaults back to rotations.
We limit our simulations to containing only two insecticides as this is the current number in mixture IRS and LLINs. We therefore F I G U R E 3 Example Simulations of the Sequence, Rotation and Mixture Strategies. Sequence strategy: Insecticide 1 (red) is deployed continuously until it reaches the withdrawal threshold and is replaced by the Insecticide 2 (blue) at the next deployment opportunity. Insecticide 2 is deployed continuously until it too reaches the withdrawal threshold. When Insecticide 1 is not deployed, fitness costs reduce resistance. Once the resistance has reached the return threshold (at 300 generations), it becomes re-available for deployment. The sequence simulation terminates at 420 generations as neither insecticide is available for deployment. Rotation Strategy: The deployed insecticide deployed is changed at each deployment interval. The simulation was terminated at 500 generations, with both insecticides still available for deployment as neither had yet reached the withdrawal threshold. Mixture Strategy: Insecticide 1 and 2 are deployed together as a single formulation, in a de facto sequence. The simulation was terminated at 500 generations, with both insecticides still available for deployment as neither had yet reached the withdrawal threshold. These simulations are presented only to explain how each IRM strategy works. The parameter inputs were: deployment interval = 10 generations, cross resistance = 0, start resistance = 0, heritability = 0.2028686, male insecticide exposure = 0.2500806, female insecticide exposure = 0.8645515, fitness cost = 0.1730878, intervention coverage = 0.6885651, dispersal = 0.8844445. look at purely evaluating the differences between strategies, and not the complexities associated with allowing for additional insecticides and extending the armory of strategies. We note that the research, development, and economic complexities with creating a mixture formulation may not be equivalent to creating two separate single insecticide formulations.

| Simulations to directly compare strategies
This set of simulations assumes each insecticide in the parameter combination has the same properties (i.e., heritability, fitness cost and starting resistance are the same for both insecticides). This allows for better comparison of the effect of the IRM strategy itself as the efficacy is not dominated by a single insecticide to which resistance evolves very slowly or inhibited by an insecticide to which resistance evolves very rapidly (this assumption is relaxed in the next section).
The parameter space is described on Table 3 and was sampled using Latin hyperspace sampling (Carnell, 2020). A total of 5000 randomly generated parameter inputs were used. Cross resistance/ selection values were used from −0.5 to 0.5 at 0.1 intervals, giving 11 different values. Simulations were started assuming both insecticides were novel (z = 0, in both intervention site and refugia) or both insecticides had previously been used and there was some resistance present (z = 50, corresponding to ~5% bioassay survival).
Therefore, a total of 220,000 (5000 random parameter values × 11 cross resistance values × two starting resistance values × two deployment intervals) unique parameter inputs were available for each IRM strategy. The same 220,000 parameter inputs were used for each IRM strategy to allow for direct comparisons between strategies. Three IRM strategies were evaluated so a total of 660,000 (220,000 × 3) simulations were conducted.

| Simulations with unique insecticide properties
We ran a set of simulations which allowed for each insecticide in the armory to have different, unique properties. In these simulations the two insecticides were allowed to have different heritability, fitness cost and starting resistance (PRS). For the starting PRS, this was sampled within a uniform distribution between 0 and 80 (ensuring the PRS was below 100 as this is our criteria for a failed insecticide). For each insecticide, the starting PRS was the same in both the intervention site and the refugia. Cross resistance was included between the insecticides and was sampled within a uniform distribution between −0.5 and 0.5. Latin hyperspace sampling was used again (Carnell, 2020

| Outcome measures
One challenge with any study of different strategies/interventions/ policies is defining how best to evaluate them (see discussions in Madgwick & Kanitz, 2022a). For evaluating IRM we consider that there are three potential main outcomes of interest. First is the "strategy lifespan", second is the average level of resistance and third is the highest level of resistance achieved. We consider the "strategy lifespan" to be the primary outcome.
Strategy Lifespan: The total duration of the simulation in generations. The simulation stops either when (i). neither insecticide is available for deployment (because both are above the return threshold) or (ii). when the simulation has run for 500 generations. The reason for capping at 500 generations is if there are no obvious difference between deployment strategies at this point, which represents a ~50-year time horizon, then the strategies are equivalent over any notional policy timeframe. When one strategy has a longer strategy lifespan than another for the same parameter inputs, we define that strategy as having won and the other strategy to have lost. If two (or more) strategies have the same strategy lifespans, we say they have drawn. We presume the longer an IRM strategy lasts, the better the IRM strategy. We report the number of wins, losses and draws as percentages of the totals. We also note the size of the win in terms of percentage difference in lifespan on the insecticides.
We expect a strategy to have a strategy lifespan of at least 10% longer to justify using a potential more logistically complex or economically expensive IRM strategy. We therefore define any win where the strategy lifespan increased by 10% or more to another as an "operationally relevant win". This is especially the case where the benefit of one IRM strategy over another may not be seen for over 30 years, which is beyond the timeframe used for operational planning. This threshold was used as a soft cut-off as used similarly by Madgwick and Kanitz (2022b).
If the strategy lifespan is ≥10% longer than another strategy (under the same conditions) we describe that strategy as having an "operationally relevant win". This helps to identify simulations where there would likely be a benefit in choosing one strategy over another.
When two strategies have equal strategy lifespans (for the same parameter inputs), we say they have drawn. Draws can occur if both simulations run out of insecticides (i.e., all insecticides have bioassays survival <10%) at the same timepoint, or a pair of simulations both reach the 500-generation maximum. Where the sets of simulations draw, the secondary outcomes can be compared. Secondary outcomes are compared only in simulations with equal strategy lifespans.
The secondary outcomes are: 1. The mean bioassay survival to the currently deployed insecticide(s) during the simulation. We report this rather than the corresponding field survival (Equation 1b) as bioassay survival is measured in operational settings to inform decisions. We assume insecticide-based interventions are more effective when bioassay survival is lower as the killing effect on the population size and age structure would be greater. Therefore, a lower mean bioassay survival is preferrable.
2. The peak bioassay survival reached during the simulation by any insecticide at any timepoint. As bioassay survival increases, we may expect an increase in the potential for compensatory genes to evolve to counteract the effect of fitness costs, and therefore a lower peak bioassay survival is preferable.

| Statistical analysis
All statistical analysis and data visualization was conducted in R version 4.0.3 (R Core Team, 2020). The following R packages were used: epiR (Stevenson et al., 2020) for partial rank correlation, mgcv (Wood, 2011) for generalized additive models, MASS (Venables & Ripley, 2002) for the negative binomial GLM, and randomForest (Liaw & Wiener, 2002) for the random forest models. Data visualization used ggplot2 (Wickham, 2009). We separate our analysis into three parts.
• First, we directly compare strategies deploying a single insecticide at a time (sequences and rotations) against one another.
• Second, we compare deploying insecticides singularly (sequences and rotations) against deploying the same insecticides together in a mixture.
• Thirdly, model sensitivity analysis was conducted using partial rank correlation, generalized linear modeling, and random forest models to identify parameters driving IR and offering predictive value in informing the IRM decision-making process.

| Partial rank correlation
Partial rank correlation was used to assess the degree of correlation between the randomly generated parameter input values and the strategy lifespans of the IRM strategies where the insecticides were given equal properties. This identified the main factors driving the evolution of IR and was conducted separately for each IRM strategy, each level of starting resistance, and each cross-selection value. Two-sided partial rank correlation was used to determine the direction of the correlation.

| Generalized linear modeling and generalized additive modeling
Simulations which ran to completion (i.e., terminated at 500 generations) were excluded from the analysis as these simulations were artificially terminated and their inclusion could lead to an underestimation of effect sizes. A total of 86,566 simulations were included after these runs were excluded. Initial exploration with a Poisson GLM, indicated the data were over-dispersed (dispersion = 5.937677, p < 2.2e−16). A negative binomial generalized linear model was therefore fitted. The model parameters of heritability, female insecticide exposure, male insecticide exposure, intervention coverage, fitness cost and dispersal were used as predictors. The IRM strategy was input as a factor. The starting resistance was converted to a factor, with 0 being defined as a "novel insecticide" and 50 being a "preused insecticide", in this set of simulations insecticides only started at either 0 or 50 PRS. The deployment interval was input as a factor.
Generalized additive models were conducted to assess for nonlinear relationships. Where notable nonlinear relationships were found (Supplement S4, Figure S4), splines were included in the model.

| Random forest models
Random forest models can then be used to identify which individual parameters provide information which gives the best predictive accuracy. Eight random forest models were fit to predict the operational outcome. Such that the models would be predicting the optimal IRM strategy using the following comparisons. This model fitting process was repeated but restricted to include only simulations where the intervention coverage was ≥0.5. when coverage is >0.5, the no operational win outcome was no longer the dominant outcome (Supplement S4, Figure S5).
The random forest models were fit against a random sample of 70% of the respective simulations (a training dataset) and was used to predict the operational outcome for the remaining 30% of the samples to estimate the accuracy of the model in predicting the correct outcome. The variable importance of each of the parameters from each model is then reported and corresponds to how much the prediction error is affected by removing the parameter (Liaw & Wiener, 2002). Therefore, parameters with greater importance correspond with higher model accuracy, and are likely to be more beneficial to measure in the field as a basis for IRM decision-making.

| Comparing sequences versus rotations
When comparing the strategy lifespans at the "global level" (i.e., crude averages, not accounting for any other parameters), then the sequence and rotation strategies appear to perform broadly equally, with a similar number of wins and operational wins for each strategy (Table 4). However, when accounting for the starting resistance, deployment interval and the degree of cross resistance more distinct patterns can be observed. When the cross resistance between the two insecticides is positive, sequences become the preferable strategy to maximize the strategy lifespan. If cross resistance is zero or negative (i.e., resistance to one enhances susceptibility to the other) rotations become the favored strategy, unless the insecticides have pre-existing resistance. In this case there is little difference with a 10-generation deployment interval, but sequences are favored with a deployment interval of thirty generations (Figure 4).

TA B L E 4 Comparing sequences and rotations.
Outcome: Considers the absolute length of the simulations Note: Sequences and rotations were each run for the same set of input variables (Table note 1 below) and deployment intervals, for a total of 220,000 sets of input variables. The winning strategy was determined as described in Table note 2 below. Variables include for example, heritability, male exposure, female exposure, starting resistances, cross-resistance values etc. In the first "absolute criterion" the winning strategy was the one with the longest time before resistance had reached the 10% withdrawal threshold for all the insecticides in the armory. A draw was recorded if both lasted the same time, or both lasted until the end of the simulation (500 generations, ~50 years). The "operational criteria" works in an analogous manner, but a "win" only occurs if a strategy lasts >10% longer than the other. The rationale behind this "operational "criterion is that if the difference is less than 10% then the choice of "best" strategy will likely depend on operational or logistic considerations.

F I G U R E 4
Frequency distributions comparing the strategy lifespans between sequences and rotations for insecticides with identical properties. Values above zero indicate rotations (red) are the favored IRM strategy, values below zero indicates sequences (blue) were the better performing strategy. Draws (identical strategy lifespans, where the difference was 0 generations) were excluded from the plots. Left plot is the overall global distribution without any stratification. The four plots on the right are stratified by the deployment interval (top-bottom: 10 or 30 generations) and starting resistance (left-right; novel (z = 0, 0% bioassay survival) or pre-used (z = 50, ~5% bioassay survival)). Each row in these plots is the frequency distribution of the percentage difference in strategy lifespan for depending on the amount of cross-resistance. For the cross-selection plots, the numbers on the right-hand side of each plot indicate how many comparisons are plotted in the distribution (that is, the number of comparisons where the simulations did not draw).
Due to the large number of draws in strategy lifespans (72% of comparisons are tied) it is useful to also consider the secondary bioassay survival outcomes. When sequences and rotations draw, rotations outperform sequences by having both lower peak bioassay survival and a lower mean bioassay survival to the deployed insecticide, regardless of the starting resistance, deployment interval, or degree of cross resistance ( Figure 5). This is expected because most draws occur when both strategies last for the maximum of 500 generations and sequences by model rule definition explicitly run insecticides up to the withdrawal threshold. In contrast, rotations often hold both insecticides at lower PRS over the duration of the simulations.

| Comparing sequences and rotations versus mixtures
Where it is chemically possible to combine two insecticides together into a single mixture insecticide formulation, should these insecticides be deployed as a single mixture formulation or as two separate single insecticide formulations? Mixtures were found never to perform worse when compared against sequences ( Figure 6) or rotations ( Figure 7) regardless of the starting resistance, deployment interval or degree of cross resistance between the two insecticides when both insecticides were given identical properties.

| Comparing sequences, rotations, adaptive rotations and mixtures with unique insecticides
When allowing the insecticides to have unique properties (starting resistance, heritability, and fitness cost), the rotation strategy performed badly. However, this is due to the deployment rules of the rotation model, which prevent immediate redeployment that is, if one insecticide has failed there is no longer anything available to be rotated so the simulation is terminated, meaning that the strategy can no longer be implemented. When the deployment rules are relaxed allowing rotations to be adaptive rotations (defaulting to be in sequence when necessary), adaptive rotations generally outperformed sequences (Figure 8). This is expected because the adaptive rotation strategy rotates but with facility to become a temporary sequence when required.
When the insecticides in the mixtures were given unique properties (starting resistance, heritability, and fitness cost), there are occasions when mixtures are no longer the dominant strategy, although these situations were rare when compared to the number F I G U R E 5 Difference in the peak and mean bioassay survival for drawn sequences versus rotations simulations for insecticides with identical properties. Left Panel is difference in the peak bioassay survival, where values above zero favor rotations and values below zero favor sequences. Right panel is difference in the mean bioassay mortality to the deployed insecticide, where values above zero favor rotations and values below zero favor sequences A total of 158,337 simulation pairs (same parameter inputs) which drew are included. In summary, rotations never appear to perform worse (when evaluated on the secondary outcomes) in the draws as they maintain the resistance at a lower level.

F I G U R E 6
Frequency distributions comparing the strategy lifespans between sequences and mixtures for insecticides with identical properties. Percentage differences above zero indicate mixtures were the favored strategy, values below zero favor sequences. Draws (identical strategy lifespans, difference = 0) were excluded from the plots. The left plot is the overall global distribution without any stratification. The four plots on the right are stratified by the deployment interval (top-bottom) and starting resistance (novel: z = 0, preused: z = 50; left-right). Each row in these plots is the frequency distribution of the percentage difference in strategy lifespan for whether cross-resistance was positive, negative, or not included. Difference values above 0 (green) favor mixtures. For the cross-selection plots, the numbers on the right-hand side of each plot indicate how many comparisons were in each distribution (that is, the number of comparisons where the simulations did not draw). In summary, mixtures never appear to be inferior to sequences when both insecticides have equivalent properties.
F I G U R E 7 Frequency distributions comparing the strategy lifespans between rotations and mixtures for insecticides with identical properties. This figure has exatly the same structure and interperetion as Figure 6, it simply compares rotations (rather than sequences) agianst mixtures.
of mixture wins (Figure 8). In these situations, the insecticides were very dissimilar from one another in terms of their starting PRS and heritability (Supplement S4, Figure S2). This indicates that only one insecticide in the mixture failed, and the other insecticide would remain available to be deployed singularly as a sequence.

| Identifying parameters driving IR to inform decision-making
There are situations where while on average two or more IRM strategies may perform equally well (Tables 4 and 5

| Random forest models
A total of eight random forest models were fitted on the operational outcome predicted by the model inputs, with the aim to identify parameters useful for predicting which strategy to use. Random forest models allowed for the identification of parameters which offer predictive power in determining which IRM strategy to use. Parameter importance is identified by the mean decrease in accuracy. Where higher values indicate these parameters are more important in the model, and therefore more important in terms of providing information on which to make informed IRM decisions. Figure 9 shows the parameter importance when choosing between sequences and rotations. Figure 9a being the parameter importance when both insecticides had equal properties, and highlights intervention coverage and the degree of cross resistance between the insecticides as important predictors. Figure 9b shows the parameter importance when both insecticides were allowed to have unique properties, here intervention coverage and cross resistance were also found to be the predictors with the greatest importance. Figure 10 shows the parameter importance when choosing between sequences, rotations, and mixtures. Figure 10a is for when both insecticides had the same properties and highlights the importance of intervention coverage followed by heritability. Figure 10b is when the insecticides had unique properties intervention coverage became more important (relative to the other predictors), with cross resistance the next most important predictor.
Intervention coverage being an important predictor is intuitive because it dictates the amount of insecticide selection. However, its use as a predictor lies more in deciding whether an IRM strategy is needed at all, because at lower intervention coverages, the "no operational win category" is dominant because all strategies are likely to run to the maximum of 500 generations (Supplement S4, Figure S5). Therefore, a second set of random forest models was fit, but including only simulations where the intervention coverage was greater than or equal to 0.5, to better identify parameters which help to choose between IRM strategies.

| DISCUSS ION
We modeled polygenic IR as a quantitative trait by devising a Polygenic Resistance Score that can be measured as WHO cylinder Note: The methods and winning criteria are as described on Table 4. The difference between this Table and Table 4 is that mixtures can beat one strategy (e.g., sequences), while losing to the other (e.g., rotations; this example would be classed as a "Sequences Lose).

TA B L E 5
Comparing mixtures against sequences and rotations.

F I G U R E 9
Parameter importance from random forest model for choosing between sequences and rotations. Panel a (blue bars) is the parameter importance for when both insecticides have equivalent properties and it has a predictive accuracy of 87.65%. Panel b (green bars) is the parameter importance when the insecticides are given unique properties and has a predictive accuracy of 90.98%. The x axis shows the decrease in predictive accuracy when the parameter is removed from the analysis so larger values indicate greater predictive importance. Parameters are ordered by importance.
bioassay survival. Previous models have applied a quantitative genetics framework to IR traits (Gardner et al., 1998;Haridas & Tenhumberg, 2018), but did not investigate the implications for resistance management strategies utilizing multiple insecticides. Gardner et al. (1998)  generations per year) and the exposure scaling factor. There are of course important caveats of the model to address to provide context to our results.

| Model caveats
The first caveats concern the rules of insecticide deployment and withdrawal used. We highlight our simulations were terminated under idealistic deployment conditions when mosquito survival is still relatively low (i.e., <10% in bioassays). Setting the withdrawal threshold higher (e.g., 20%, 30% or 50% bioassay survival), would inevitably lead to more strategies drawing (as they would run to 500 generation completion), and requiring comparison based only on the secondary outcomes. Or setting a higher withdrawal threshold would require re-calibrating the simulations such that the exposure scaling factor is increased also to maintain an average 10-year insecticide lifespan under continuous deployment.
Insecticide deployments in the models are based on clearly de- reach the 10% bioassay threshold in ten-generation deployment interval five generations after deployment, then the simulation terminates five generations later. However, if the deployment interval was 30 generations, the simulation would continue for a further 25 generations before terminating, artificially inflating the apparent strategy lifespan. We feel this reflects the operational reality as LLINs are designed for 3-year deployments (~30 generations) and IRS for 1- year deployments (~10 generations) and it is unlikely that detection of resistance would result in replacement of nets after 1.5 years or IRS after 6 months. Insecticides are also deployed "perfectly" in our models, and we do not account for economic constraints dictating the continued use of cheaper insecticides. This is especially important when considering mixtures, as these formulations are expected to be more expensive and there will be an economic temptation to reduce the concentration of each insecticide in the mixture. Neither do we allow for any logistical constraints with regard to insecticide deployment. Insecticide deployments always occur instantaneously at the set time (without the delays often noted in practice).
Insecticide choices are also made "perfectly", the resistance in the field is tracked without error and information is available to make instantaneous, correct decisions regarding insecticide deployment.
There is likely to be a considerable lag between insecticide susceptible testing, insecticide procurement and insecticide deployment resulting from inadequate surveillance systems and the release of procurement funds (Chanda et al., 2015). In operational reality, insecticides are highly unlikely to be deployed in the idealistic ways the rules the model demands. It is also unlikely that a single IRM strategy would remain fixed for eternity, with an amalgamation of sequences, rotations and mixtures potentially being used in a single intervention site. We could only realistically investigate the IRM strategies under these idealized conditions of "perfect" deployment because space restraints prevented a more extended investigation of how operation factors may alter the relative merits (although the model can later be used to investigate specific examples of operational problems). We can however speculate on how operation limitations may affect the strategies. Mixtures are likely to be the most robust because no replacement decisions are required (Figure 3) supporting our conclusion that full-dose mixtures are overall the "best" IRM strategy. Sequences and rotations performed very similarly in our simulations but as noted by Curtis (1987) andHastings et al. (2022). the fact that rotations have preplanned insecticide replacement already in place may make them more operationally robust than sequences whose unpredictable replacement intervals require rapid responses.
Second are caveats concerning nonpublic health insecticide selection pressure. The model only tracks insecticides in the armory

F I G U R E 11
Parameter importance from random forest model for choosing between sequences and rotations restricted to simulations where intervention coverage ≥0.5. Random forest model was restricted to comparisons where intervention coverage ≥0.5. Panel a (blue bars) is the parameter importance for when both insecticides have equivalent properties and has a predictive accuracy of 79.09%. Panel b (green bars) is the parameter importance when the insecticides are given unique properties and has a predictive accuracy of 84.73%. The x axis shows the decrease in predictive accuracy when the parameter is removed from the analysis so larger values indicate greater predictive importance.
and assumes these insecticides are only ever deployed (and therefore encountered by mosquitoes) as part of a public heath deployment. Other insecticides, which mosquitoes may encounter such as agricultural insecticides or personal household insecticides (i.e., coils and aerosols) were not considered in the model as they are not part of the armory and we assumed they had no cross resistance with insecticides in the armory. Agricultural insecticide use is associated with increased IR in malaria vectors (Abuelmaali et al., 2013;Nkya et al., 2014;Reid & McKenzie, 2016). This further extends to the model's intervention site and refugia existing "isolated" from rest of the world. It would be possible to build these "external" selection pressures into our model (through Equations 3a and 3b) and to allow them to act even in refugia (by adding an additional term representing external selection in Equation 9a) but this is a specific operational question we could not realistically address in this study due to time constraints.

| Model interpretation
With these caveats in mind, we can say, from a purely IRM perspective the overall best strategy appears to be full-dose mixtures (Figures 6 and 7). This broadly corroborates with the results from monogenic models (e.g., Madgwick & Kanitz, 2022b;South & Hastings, 2018) for favoring full-dose mixtures. There appears to be little difference between rotations and sequences when comparing operational longevity at the "global level" (Figure 4) agreeing with monogenic modeling results (Hastings et al., 2022), although rotations are slightly more beneficial when considering mean and peak bioassay survival ( Figure 5). The reason for this would appear to be that sequences are in practice a type of rotation Hastings et al. (2022) that is, because the first insecticide in the sequence can later be re-deployed if resistance has fallen below its redeployment threshold. A blend of these two strategies (which we defined as adaptive rotations) generally was the better strategy when considering solo insecticide deployments (Figure 8).
These conclusions ignore important operational and economic factors required for the IRM implementation, for example the increased costs of a mixture LLIN/IRS, the infrastructure and supply chains required for annual rotations, and so on. However, these results highlight a vital conclusion from our study that is, the underlying assumption of a purely monogenic or purely polygenic basis of IR does not overly impact the general conclusions from mathematical models; lending further support to mixtures being the generally optimal IRM strategy. Validating IRM models with field data is challenging, probably impossible, because of the need for long field surveillance of the levels of resistance the need to measure other parameters associated with quantifying the degree of selection, and the need for replication over several IRM sites to draw robust empirical conclusions. Therefore, mathematical modeling studies using different approaches and assumptions yet yielding qualitatively similarly conclusions suggests there is some robustness to the conclusions made.
One advantage of taking a quantitative genetics approach is the ability to easily incorporate cross resistance/selection using correlated responses. One important concern with mixtures has been the role cross resistance could play (e.g Curtis, 1985).
Importantly, the presence of cross resistance between insecticides in the mixture does not appear to invalidate the conclusion that mixtures are the more beneficial IRM strategy (Figures 6 and 7), despite previous speculation cross-resistance would compromise mixtures as an IRM strategy. There are two possible explanations for this. First, because cross resistance also reduces the strategy lifespan of sequences and rotations through indirect selection occurring on the non-deployed insecticide. Second, our mixture consists of both insecticides at full-dose and this double-dose effect likely overcomes the detrimental impacts of cross resistance. In the few examples where a monogenic model has included cross resistance this pattern was also observed (Birch & Shaw, 1997;Roush, 1998;Sudo et al., 2018). This is certainly not to say we (nor indeed the authors of the previous studies) recommend mixing two insecticides where cross-resistance or cross-selection is likely, but that if these two insecticides were to be deployed in the same area mixtures would very likely be a better IRM strategy deploying them singularly as sequences or rotations. We would instead emphasize the importance of insecticides in the entire insecticide armory being chosen that have as a little cross-resistance as possible. Note, the simulations were terminated when mosquito survival was still relatively low (i.e., <10% in bioassays). At higher resistance levels mosquito survival increases, and the protective effect of each mixture component therefore decreases. Whether cross-resistance becomes more detrimental under lower insecticide doses and higher resistances needs investigation, as these are both real-world problems.
TA B L E 6 Summary of results presented in the figures. The model is calibrated such that an "average" insecticide lasts 10 years 3 Explains how insecticides are deployed and the concept of the withdrawal and return thresholds 4 Sequences and Rotations perform equally at the "global level" Rotation better than sequences negative cross-resistance Sequences slightly better than rotations with positive cross-resistance 5 When sequences and rotations draw on strategy lifespan rotations have: • Lower peak bioassay survival • Lower mean bioassay survival to deployed insecticide. Mathematical modeling allows insights into how decisions could or should be made from a policy or operational perspective, and also identifies important research and knowledge gaps through the identification of understudied yet operationally important variables. Our analysis highlights the importance of cross resistance between insecticides, intervention coverage, IR heritability and the fitness costs associated with IR as being important predictors for which IRM strategy to deploy (Figure 12). Unfortunately, the heritability of resistance is challenging to reliably measure in the laboratory (Rosenheim, 1991) and, in any case, such estimates are unlikely to reflect heritability in the field where environmental variance will be much higher. Fitness costs are also challenging, as these can affect a wide range of lifehistory traits and are difficult to measure outside unrealistic laboratory settings. There is also large variation in the design and execution of fitness studies (Freeman et al., 2021)  There is a need to manage the longer-term effectiveness of insecticides alongside the immediate epidemiological benefits of deploying insecticides and we argue that computer simulation is the most feasible way to identify best-practice in deployment of insecticides used for vector control.

ACK N OWLED G M ENTS
We thank Andy South for comments on early drafts and coding advice.

FU N D I N G I N FO R M ATI O N
NH is an MRC-DTP candidate. The funder had no role in the study design, data collection, interpretation of data, writing of paper or decision to publish.

CO N FLI C T O F I NTER E S T S TATEM ENT
The authors declare they have no competing interests.

DATA AVA I L A B I L I T Y S TAT E M E N T
Model code for running and analyzing simulations is available from https://github.com/NeilH obbs/polyres and/or the corresponding author on request.