In ian-hastings/rotations: Simulating the Evolution of Insecticide Resistance under Rotations

rotations_notes.Rmd

4/10/2017

How to make the dataframe that saves results in run_rot flexible to the number of insecticides ?

df_results <- data.frame(generation=1:max_gen, insecticide=NA, r1_active=NA, r1_refuge=NA, r2_active=NA, r2_refuge=NA, r3_active=NA, r3_refuge=NA, r4_active=NA, r4_refuge=NA, stringsAsFactors = FALSE)

Can I adapt this from array_named() in resistance ?

' n_insecticides <- 3

' # list(c()) is important below to create one list for each dim

' l_dimnames <- rep(list(c('no','lo','hi')), n_insecticides)

' names(l_dimnames) <- paste0("niche", 1:n_insecticides)

' do.call(array_named,l_dimnames)

this feels hacky, but it does work & I can put in a function

create the flexible r* columns

n_insecticides <- 3 l_ <- rep(list(NA), n_insecticides*2) #2 because active & refuge names(l_) <- c( paste0("r", 1:n_insecticides, "active"), paste0("r", 1:n_insecticides, "_refuge") ) df_r <- do.call(data.frame,l)

create inflexible columns ready to add r* columns on

max_gen <- 10 df_results <- data.frame(generation=1:max_gen, insecticide=NA)

df_results <- data.frame(df_results, df_r)

Note this can only cope with one insceticide being active at a time ? (I could put a comma delimited string in later if I did want to record multiple insecticides).

Does this fit with tidy data ?

Each variable you measure should be in one column. Each different observation of that variable should be in a different row.

Yes I think it does. Although maybe not ? Currently I record df_results and gather to df_r2 to allow plotting. To get : generation insecticide region resistance 1 1 NA r1_refuge 0.001000000 2 2 1 r1_refuge 0.001061247

So perhaps I can just save as that to start ? It has multiple rows per generation. More difficult to see, but easier to plot ?

This allows me to facet by region. BUT if I wanted to put the lines for active & region into the same subplots (which would make easier to view) I would probably want columns for active and refuge or would I ??

gen active_insecticide resistance_gene active_or_refuge resistance

I can use tidyr::separate() to get from r1_refuge to r1 & refuge in different columns.

separate(region, into=c("resistance_gene","active_or_refuge"))

I now have run_rot() working for flexible numbers of insecticides (1-5).

done : Next, get active and refuge onto the same plots.
done : allowed starting frequencies to be set

implemented set_expsoure_rot() now working dfr <- run_rot(rot_interval=50, expo_hi=0.9) dfr <- run_rot(rot_interval=50, expo_hi=0.2)

diff frequency curves for r1 & r2

dfr <- run_rot(n=2, rot_interval=50, expo_hi=c(1,0.1))

started fitness_single_locus() from run_rot()

dfr <- run_rot(n=2, rot_interval=50, eff=c(1,0.1))

9/10/17 ~ started first shiny UI : rotresist1 ~ add rot interval and generations to rotresist ~ looking good

10/10/17

~ seems that resistance evolves more slowly in later insecticides despite same params. Why ?

~ aha! the new UI helped me to see that this can be a result of the decline in resistance frequency for the later insecticides when they are not being used at the start. e.g.

run_rot( n_insecticides = 3 , max_gen = 150 , start_freqs = 0.1 , rot_interval = 50 , eff = 0.6 , dom = 0.5 , rr = 0.5 , expo_hi = 0.5 , coverage = 0.8 , migration_rate_intervention = 0.01 , cost = 0.1 )

When cost is set to 0 all 3 insecticides give the same curves :

done ~ in plots change r1 to insecticide_1 etc. done ~ in plots add the deployment symbol to legend

done ~ can I offer option to have plots on log y axis to make it clearer whats going on at lower frequencies ?

Now from the UI ticking log y axis influences the next plot, (but is isolated so allows you to compare log & non log plots).

Starting to think about sensitivity analysis. For Levick et al the inputs and ranges were hardcoded into resistance::sensiAnPaperPart()

For rotations

dominance of selection | dominance of cost | 0 - 1

Sudo manuscript rotations are always a single generation.

How do we deal with the problem that resistance frequenciesin the rotation scenario may go above 50% and this makes time-to-resistance not a good measure of how good the intervention is (because the rotation).

Ian wants it to be a fair comparison.

I'm nervous about making the rotation responsive. We both said that i unlikely to happen in the field.

Could we have number of generations at < 50% ? Not perfect either.

~ setting dominance of cost and dominance of selection to be different in fitness_single_locus().

31/10/17

UI : Warning: Error in run_rot: argument 6 matches multiple formal arguments

probably due to recent splitting of dominance

Yes, fixed now.

Increasing dominance of cost decreases resistance, increasing dominance of selection increases resistance. (both as expected)

Ians document : manuscript notes.doc

We will compare 4 strategies • RwR (rotate when resistant.. really a sequential use policy but the term RwR emphasises that can return to insecticides used earlier in the sequence if they are now effective • Periodically rotate every 10 generations (provided there replacement insecticide has resistance levels below the threshold. 10 generations is about every year • Rotate every 20 generations (~ 2 years) • Rotate every 50 generations (~ 5 years).

Threshold for a policy change: resistance allele frequency > 50% (may also want to investigate a lower level e.g. 25% and/or a fitness criterion)

End of policy life: when resistance levels to all insecticides exceed the threshold for policy change.

Sensitivity analysis to test which of the policies is best will explore the parameter space defined on table 1. However, there are a few test cases we will examine first to test our understanding and illustrate the basic principles.

Test 1. No refugia (i.e. Coverage=100%) or fitness costs. Would expect all policies to be equivalent.

Test 2. Dominance of resistance and fitness costs both set to 0.5. I suspect one of the advantages of rotations (and the presence of refugia) is that they maintain lower resistance allele frequencies which slows the spread of recessive resistance alleles. I suspect/guess/intuit that if we make them semi-dominate these effects should disappear and all policies should be the same. I’m probably wrong but it’s a nice test…..

Things still to consider RAF can go below the original defined level due to natural selection…de we want to stop this e.g. by saying this minimum value is mutation/selection

What should we do about planned rotations where RAF>50% for several generation i.e. is an ineffective deployment. This gives it an ‘unfair’ advantage over RwR

We still need to set a trap so can’t rotate on the basis of a few generations e.g. every 2 gens indefinitely as selection brings RAF>0.5 and natural selection brings RAF<50%

Notes: List of research questions (draft keep adding)

• What are the relative impacts of fitness costs and size of refugia. Presumably we find lifespan of policy then do PRC (partial rank correlation) to quantify their impact.

• Routine rotation may be better than a “rotate-when-resistant” (RwR) policy if resistance is recessive i.e. rapid rotation ensures resistance allele frequency in treated populations remains so low that resistant homozygotes are rare.

• Can contrast to mixtures. Would have to use the full model for mixtures and this one for rotation. Unfortunately, we can’t have a separate refugia in the full model at present but will build this when simulating macro-mosaics. Once that is done we can use the full model for rotations so don’t need the simple model described here, although the algebra above may be faster than running the full model.

We are not necessarily recommending rotate when resistant (RwR) policy, but it serves a good baseline to compare reactive policies (such as RwR) with pro-active polices such as rotations.

RwR may be good if there is an exceptionally long lived insecticide in the armoury as will find and use it. If that long-lived insecticide is part of a rotation sequence, you may never find it (until it has outlasted all the others).

When RwR you know how many insecticides you have left. In a strict rotation, they may all start to exhibit resistance at roughly the same time with little fore-warning.

~ looking at fixing migration rate to be between 0 & 1, where 1 is random movement.

~~ currently suggestion is that it is set to (0.1-0.9)*(1- Coverage)

So I want to remove dependence on the 1-coverage bit, by putting that into the code.

migration_rate_intervention = 0.01,

migration rate into and out-of the treated area. It is the proportion of the treated population that migrates. We assume that immigration=emigration.

First line of run_rot()

migration_rate_refugia = migration_rate_intervention * coverage/(1-coverage)

This is done for all insecticides

 fem_intervention=
   (1-migration_rate_intervention)*RAF[temp_int, 'f', 'intervention', gen]+
   migration_rate_intervention*RAF[temp_int, 'f', 'refugia', gen]

 male_intervention=
   (1-migration_rate_intervention)*RAF[temp_int, 'm', 'intervention', gen]+
   migration_rate_intervention*RAF[temp_int, 'm', 'refugia', gen]

 fem_refugia=
   (1-migration_rate_refugia)*RAF[temp_int, 'f', 'refugia', gen]+
   migration_rate_refugia*RAF[temp_int, 'f', 'intervention', gen]

 male_refugia=
   (1-migration_rate_refugia)*RAF[temp_int, 'm', 'refugia', gen]+
   migration_rate_refugia*RAF[temp_int, 'm', 'intervention', gen]

RAF[temp_int, 'f', 'intervention', gen]=fem_intervention
RAF[temp_int, 'm', 'intervention', gen]=male_intervention 
RAF[temp_int, 'f', 'refugia', gen]=fem_refugia
RAF[temp_int, 'm', 'refugia', gen]=male_refugia

Perhaps I can put this into a little spreadsheet to test ? Or into a little test function. e.g. that when migration set to 1 mixing is random ? Just for 1 gender.

coverage

migration_refugia = migration_intervention * coverage/(1-coverage)

freq_new <- (1-migration)*

Nearly there, testing Ians migration :

Ah yes 1-c is max migration, it gives equal frequencies in one generation

migration_test(migration=0.4, startfreqs=c(0.8,0.2), coverage = 0.6) gen site 1 2 3 4 5 6 7 8 9 10 intervention 0.8 0.56 0.56 0.56 0.56 0.56 0.56 0.56 0.56 0.56 refugia 0.2 0.56 0.56 0.56 0.56 0.56 0.56 0.56 0.56 0.56

done ~ fix migration rate to be between 0 & 1, where 1 is random movement.

~ allowed plots not to show refuge

high cost, no migration

run_rot( n_insecticides = 3 , max_gen = 150 , start_freqs = 0.01 , rot_interval = 10 , eff = 1 , dom_sel = 0.5 , dom_cos = 0.5 , rr = 0.5 , expo_hi = 0.5 , coverage = 0.8 , migration = 0 , cost = 0.1 , logy = FALSE )

no cost, no migration

no cost and rotate when resistant, no migration

run_rot( n_insecticides = 3 , max_gen = 150 , start_freqs = 0.01 , rot_interval = 0 , eff = 1 , dom_sel = 0.5 , dom_cos = 0.5 , rr = 0.5 , expo_hi = 0.5 , coverage = 0.8 , migration = 0 , cost = 0 , logy = FALSE )

no cost, some migration

run_rot( n_insecticides = 3 , max_gen = 150 , start_freqs = 0.01 , rot_interval = 10 , eff = 1 , dom_sel = 0.5 , dom_cos = 0.5 , rr = 0.5 , expo_hi = 0.5 , coverage = 0.8 , migration = 0.05 , cost = 0 , logy = FALSE )

done ~ improve plot legend to show when insecticide in use

20/11/17

~ looking at refactoring the migration code in run_rot()

Can do this to get the allele frequencies for generation 2 (by )

RAF[,,,2]

~ currently it does 4 very similar calcs for m|f, refuge|active.

~ i can replace with single, except need to be acreful not to change values before they have been used in the other calculations.

FROM : for(temp_int in 1:n_insecticides){

 fem_intervention=
   (1-migration_rate_intervention)*RAF[temp_int, 'f', 'intervention', gen]+
   migration_rate_intervention*RAF[temp_int, 'f', 'refugia', gen]

 male_intervention=
   (1-migration_rate_intervention)*RAF[temp_int, 'm', 'intervention', gen]+
   migration_rate_intervention*RAF[temp_int, 'm', 'refugia', gen]

 fem_refugia=
   (1-migration_rate_refugia)*RAF[temp_int, 'f', 'refugia', gen]+
   migration_rate_refugia*RAF[temp_int, 'f', 'intervention', gen]

 male_refugia=
   (1-migration_rate_refugia)*RAF[temp_int, 'm', 'refugia', gen]+
   migration_rate_refugia*RAF[temp_int, 'm', 'intervention', gen]

RAF[temp_int, 'f', 'intervention', gen]=fem_intervention
RAF[temp_int, 'm', 'intervention', gen]=male_intervention 
RAF[temp_int, 'f', 'refugia', gen]=fem_refugia
RAF[temp_int, 'm', 'refugia', gen]=male_refugia

}

TO :

mig_intervention <- (1-migration_rate_intervention)*RAF[,, 'intervention', gen]+
                       migration_rate_intervention *RAF[,, 'refugia', gen]

mig_refugia      <- (1-migration_rate_refugia)*RAF[,, 'refugia', gen]+
                       migration_rate_refugia *RAF[,, 'intervention', gen]

RAF[,, 'intervention', gen] <- mig_intervention
RAF[,, 'refugia', gen] <- mig_refugia

Which is more concise and shows more clearly what is going on.

testing multiplying RAF arrays

to allow making code more concise

RAF <- set_start_freqs(n=3, max_gen = 1, freqs=c(0.1,0.01,0.001)) RAF2 <- set_start_freqs(n=3, max_gen = 1, freqs=c(1,2,3))

RAF * RAF2 sex insecticide m f 1 0.100 0.100 2 0.020 0.020 3 0.003 0.003

done ~ put migration code into a function started ~ add a test that checks migration for diff insecticides & mf

~ rdpeng suggests not having non standard evaluation in package code : If you have any non-standard evaluation in your package code (which you’ll notice because of the “no visible bindings” warnings you’ll get when you check the package), go through and change any instances to use standard evaluation alternatives. This change prevents these warnings when you check your package and will also ensure that the functions behave like you expect them to when they are run by other users in their own R sessions. ~ https://bookdown.org/rdpeng/RProgDA/non-standard-evaluation.html

programing with dplyr

http://dplyr.tidyverse.org/articles/programming.html

21/11/17 ~ is it worth creating a get_freq() or getraf() function similar to in rtsetse but not as complex ? ~ pass RAF ~ get_freq(insecticide=1, site='intervention', sex='f', gen=2 ) ~ assume all for any not specified ~ have an as_df arg to return as a dataframe

ARG! got migration seeming like it's working

BUT now have come across non-standard-evaluation error in run_rot()

Fix this : Error in filter(., active_or_refuge == "active") : object 'active_or_refuge' not found

Draft that is sposed to help explain : http://rpubs.com/lionel-/programming-draft and posts in here : https://community.rstudio.com/t/should-tidyeval-be-abandoned/2238/37

I could remove the generations under 50 bit for now to try to get run_rot() working again.

removed add_gens_under50 and now run_rot() does run.

Ians text removed from the insecticide switching portion of code

#NOTE. At the moment if all but one of the insecticides are over the threshold frequency, the one under the threshold #will be repeatedly deployed. For example, if there are 3 insecticides and rotation is every 10 generations. #if #3 is being deployed and due to rotate out at generetion 100: if only #3 is under the threshold #it will not rotate out, but will continue to be used util the next scheduled rotation at which point #(a) it will continue to be use if it is the only one below the threshold #(b) it will be replaced if one of the other insecticides is now below the threshold due to fitness effects or migration #(c) the simulation will terminate if all insecticides exceed the threshold

22/11/17 ~ added sentence about interval 0 RwR to the UI

~ removed this commented out alternative to df_results #experimental # df_res_active <- data.frame(generation=rep(1:max_gen,n_insecticides), # insecticide=NA, # region=NA, # resistance=NA, stringsAsFactors = FALSE) # df_res_refuge <- df_res_active

~ take the insecticide switching out into its own function(s) insecticide_check() insecticide_switch()

~ rename rot_count to gens_this_insecticide

~ add in an option to limit switch back to previous insecticides ~ added something to put a lower limit on length of rotations for rwr

~ added option to stop resistance frequencies going below starting value (the thinking being that they could be maintained at this value by selection). It happens in the simulation when there is a cost and the insecticide is not in use.

~ deployed to shinyapps : working well now https://andysouth.shinyapps.io/rotresist1/

23/11/17 ~ sent draft model-data gap workshop invite to IVCC

~ tweaked UI input defaults to show a comparison of fixed rotation versus sequence (rwr)

start protocol for sensitivity analysis.

Currently run_rot() outputs : 'data.frame': 1600 obs. of 5 variables: $ generation : int 1 2 3 4 5 6 7 8 9 10 ... $ insecticide : num NA 1 1 1 1 1 1 1 1 1 ... $ resist_gene : chr "insecticide1" "insecticide1" "insecticide1" "insecticide1" ... $ active_or_refuge: chr "active" "active" "active" "active" ... $ resistance : num 0.001 0.00141 0.002 0.00283 0.00399 ...

insecticide : indicates which is in use resist_gene : which gene does this frequency relate to

from run_rot() now commented out

# calculate number generations under 50% resistance

df_res2 <- df_res2 %>% dplyr::filter(active_or_refuge=='active') %>% group_by(resist_gene) %>% summarise(gens_under50 = sum(resistance < 0.5, na.rm=TRUE)) %>% ungroup() %>% left_join(df_res2, by='resist_gene')

tst <- run_rot()

tst %>% dplyr::filter(active_or_refuge=='active') %>% group_by(resist_gene) %>% summarise(gens_under50 = sum(resistance < 0.5, na.rm=TRUE))

This gives generations under 50% resistance for each insecticide in the active area

1 insecticide1 134 2 insecticide2 144 3 insecticide3 154 4 insecticide4 164

I could then calculate the mean gens_under50 across all the insecticides.

tst %>% dplyr::filter(active_or_refuge=='active') %>% group_by(resist_gene) %>% summarise(gens_under50 = sum(resistance < 0.5, na.rm=TRUE)) %>% summarise(mean_gens_under50 = mean(gens_under50))

mean_gens_under50 1 149

Look back at the resistance sensitivity code and see how much of it I can use.

resistance files to do runs

sensiAnPaper1All.Rmd : code to run whole sensitivity analysis & produce plots
sensiAnPaperPart.r : sets up param values and runs for one treatment (e.g. mixture, insecticide1)
setInputOneScenario.r: creates input object for each run
runModel2.r : runs model for the passed scenarios

I would like to modify setInputOneScenario.r so I'm not referencing inputs by their location in the matrix which seems dangerous, but I might not be able to get around ?

Can I store the inputs in a dataframe with a name column, or with rownames ? Or probably it should be a dataframe with a column name for each input, this will make it easy to add multiple runs together.

resistance::runModel2() just accepts an input matrix (i.e. not individually specified params) you need to call setInputOneScenario() first. input <- setInputOneScenario() listOut <- runModel2(input)

rotations::run_rot() does currently accept individually specified params

I want to make rotations reasonably similar to resistance but make improvements where I can.

I can possibly 1. create a set_input_one_scenario() function that returns a dataframe (or list ?) of the inputs. 1. continue to allow run_rot() to accept individual inputs with an optional list

maybe a diff way of doing too using ... operator, if I have as a list can I pass directly without having to write them all out ?

in resistance this function passes all it's arguments onto another.

resistSimple <- function(...) { input <- setInputOneScenario(...)

can I create a list and then pass that to run_rot() (yes using do.call) and also save as a dataframe ? This use of do.call given in Hadleys advanced R.

linputs <- list(max_gen = 50, n_insecticides = 2) do.call(run_rot, linputs) dfinputs <- as.data.frame(linputs)

This is good in that I don't have to specify all the inputs. Also easy afterwards to rbind dfinputs together.

Can I make the rotations workflow with fewer files than the resistance one ? : rotations workflow : 1. sensi_an_rotation1.Rmd : run whole sensi analysis including setting input ranges, & analysing outputs 1. run_rot.r : a single model run

resistance workflow 1. sensiAnPaper1All.Rmd : code to run whole sensitivity analysis & produce plots 1. sensiAnPaperPart.r : sets up param values and runs for one treatment (e.g. mixture, insecticide1) 1. setInputOneScenario.r: creates input object for each run 1. runModel2.r : runs model for the passed scenarios

But may make sense to keep closer to the resistance workflow.

27/11/17

~ talked to Ian fine to send out invite for IVCC workshop now ~ ask Hilary & Martin if they have any other staff that can come ~ do a doodle poll for dates

~ checked with Ian, fine that insecticide switch criterion is female only

Q for Ian :

Is it mean generations under 50 for all insecticides that we want ?

discussion

Probably just for the currently active insecticide. Might be called total effectively controlled generations.

but then how do we compare multiple insecticides.

We might also be interested in the other insecticides recovering.

So we might want more than one output.

Output1 : Total well controlled generations. For each generation count whether the active insecticide is below the resistance threshold of 50%.

Output2 : (may not be as useful) Mean number of generations below resistance threshold across all insecticides at all times.

Scenarios.

Ian not that keen on my minimum timespan for rotate-when-resistant, because that potentially forces continued use of an insecticide that is ineffective. He would rather that something that stops switching back to an insecticide which has been used less than e.g. 5 generations ago. But I asked whether testing for resistance is likely to be tht frequent ? If testing is only once per year than changes can perhaps only be made every 10 generations.

We agreed that we could come back to this later, because it only really becomes relevant when there are fairly high costs that being resistance down when insecticides are not in use.

gens_dep_under50 : is the tot_gens_dep_under50 : summed across all insecticides

Good meeting with Rosemary about spray data.

28/11/17 Planning IVCC workshop. potential dates : 8,9,11,12 and 29th January are best.

8th is a monday I'd like to miss.

2 hour or half day ? I'm inclined to 2 hours. More likely to get people.

Sent invite.

started a word doc with outline

Decisions to make before : ~ Focus on mix v sequence, or rotations, or both

29/11/17 From WHO World Malaria report 2017 Q&A. About insecticide resistance & malria burden : http://www.who.int/malaria/media/world-malaria-report-2017-qa/en/ A large multicountry evaluation coordinated by WHO between 2011 and 2016 showed that insecticides continue to be an effective tool for malaria prevention, even in areas where mosquitoes have developed resistance to pyrethroids. Further, there is no clear correlation between insecticide resistance and trends in malaria burden: some countries with resistance to pyrethroids have shown reductions in disease burden; others with less resistance have shown an increase. Ultimately, we need more data on the public health impact of insecticide resistance and further investments in this area.

Also IRS coverage gone down probably because switching to more expensive insecticides :

The proportion of the population at risk protected by indoor residual spraying (IRS), an intervention that involves spraying the inside walls of dwellings with insecticides, declined globally from a peak of 5.8% in 2010 to 2.9% in 2016, with decreases seen in all WHO regions. The number of people protected in 2010 was 180 million globally, reducing to about 100 million in 2016. While there is yet no clear evidence of why this is happening, the decreases in coverage with IRS insecticides are occurring as countries change or rotate the chemicals they are using to try to prevent the spread of mosquitoes resistant to pyrethroids. Because alternative insecticides are more expensive than pyrethroids, increased funding for IRS is needed to maintain high coverage levels.

8/1/2018 assessing where I am with rotations

sensi_an_rotations1.Rmd is starting to do what I want.

When I set coverage to 1 and cost to 0 I get this :

run_rot(coverage=1) Error in if (RAF[insecticide, sex, site, gen] < start_freqs[insecticide]) RAF[insecticide, : missing value where TRUE/FALSE needed

I suspect because when coverage=1 site='refugia' is not populated.

BEWARE I've added a few changes with major implications in to try to fix this ...

Might want to just put a condition in to not run migration when coverage=1

# migration between refugia and intervention site
RAF[,,,gen] <- rot_migrate(RAF[,,,gen], migration=migration, coverage=coverage)

9/1/2018

getting there with model working with no refugia this works run_rot(coverage=1)

But knitting sensi_an_rotations1.Rmd gives error in RAF1gen[] incorrect number of dimension calls from insecticide_check() and I suspect because the site dimension has been dropped because it is just length 1

added drop=FALSE in insecticide_check() but still seem to get error

may need to add drop=FALSE at other places in or maybe the build hasn't worked yet

I think working now.

12/1/2018

run_rot(coverage=1, rot_interval=0)

Error in RAF1gen[current_insecticide, "f", "intervention"] : incorrect number of dimensions Called from: insecticide_check(RAF1gen = RAF[, , , gen, drop = FALSE]

19/1/2018

model is now working with a coverage of 1

From GPIRM : Idea that cost of resistance likely to decline when frequency of resistance gene has been high. (We don't have this in our models). Insecticide resistance management strategies must, however, be implemented before the resistance gene becomes common and stable in the population; otherwise, the resistant gene will not recede even if use of the insecticide causing selection pressure is discontinued. As demonstrated by a study of blowflies by McKenzie and Whitten in 1982 (3), fitness cost is not an intrinsic property of the gene. Therefore, if that gene is allowed sufficient time to become common in a population, the rest of the genome will adapt to incorporate it without a significant fitness cost. At this point, even if the selection pressure is removed, the resistance gene will remain in the population.

But I looked at the 1982 paper and not sure it supports that ? Seems to me that McKenzie1982 shows that relative fitness of the susceptible genotypes increase over time as would be expected by a decline in the concentration of the active ingredient ? (and similar to our WoS paper)

Example of sequential use and cost from GPIRM Example in Colombia: resistance genes disappeared from the vector population when selection pressure was removed. A report submitted by WHO AMRO regional entomologists showed that, in 2005–2006, resistance to pyrethroids and DDT was identified in An. darlingi. A decision was quickly made to change to fenitrothion, an organophosphate with a different mode of action, for IRS. Rapid implementation of this alternative, which thereby removed the selection pressure, reduced the frequency of resistance. In 2010, susceptibility tests showed that the frequency of resistance genes in the vector population had dropped below the level of detection, and pyrethroids were once again introduced into the IRS programme, albeit on a more limited scale.

Working on rotations_story.Rmd to create visuals about rotations we can use in upcoming talks.

22/1/2018

Why do rotations last longer than sequences ? (even when coverage=1 & cost=0)

I think it's just because the rotation interval continues to its end.

e.g. added an option to plot a 2nd scenario to rot_plot_resistance()

dfrot0 <- run_rot(n_insecticides=2, rot_interval=0, expo_hi=0.5, eff=0.5, cost=0, coverage=1) dfrot15 <- run_rot(n_insecticides=2, rot_interval=15, expo_hi=0.5, eff=0.5, cost=0, coverage=1) rot_plot_resistance(dfrot0, df_resanother=dfrot15, plot_refuge=FALSE)

So maybe we should have assessment of the resistance frequency in the rotation strategy too ? I could add it as a flag.

done, now works for simple scenario

BUT rotations still lasting a few generations longer in later scenarios ? dfseq <- run_rot(n_insecticides=2, rot_interval=0, expo_hi=0.5, eff=0.5, cost=0.01, coverage=1) dfrot <- run_rot(n_insecticides=2, rot_interval=15, expo_hi=0.5, eff=0.5, cost=0.01, coverage=1) rot_plot_resistance(dfrot, df_resanother=dfseq, plot_refuge=FALSE)

rotation insecticide2 192 2 active 0.46672149 insecticide2 193 2 active 0.48524932 insecticide2 194 NA active 0.50379224

sequence insecticide1 187 0.4710576 insecticide1 188 0.4895936 insecticide1 189 0.5081326

Could it be because an extra generation is slipped in at each rotation switchover ? Or might be that more decline happens in rotations

Actually I think it's due to min_rwr_interval which was set at 5 by default and in the current code is implemented for sequence but not for rotation.

hmmm but even when I set min_rwr_interval to 1 there seems to be the difference ? is there any info in the console messages ?

rotation :

generation 166, switch from insecticide 1 to 2; frequencies = 0.402193 and 0.178446 generation 181, switch from insecticide 2 to 1; frequencies = 0.402193 and 0.384221 generation 188, switch from insecticide 1 to 2; frequencies = 0.512362 and 0.393770

simulation terminating at generation 194 because all RAFs above threshold of 0.500000

frequency of resistance in females to insecticide 1 is 0.504826 frequency of resistance in females to insecticide 2 is 0.503792

sequence :

generation 92, switch from insecticide 1 to 2; frequencies = 0.510798 and 0.001000 generation 183, switch from insecticide 2 to 1; frequencies = 0.510798 and 0.397957

simulation terminating at generation 189 because all RAFs above threshold of 0.500000

frequency of resistance in females to insecticide 1 is 0.508133 frequency of resistance in females to insecticide 2 is 0.503261

Hi Ian,

I'm looking into rotations and mostly I find very little difference from sequences even with costs or refuges. In a sequence resistance builds up once and then declines, in rotations it builds up in steps and declines in between but the time-to-resistance is the same for both.

Then I refound this idea in GPIRM that costs of resistance are likely to decline when frequency of resistance gene has been high.

"Insecticide resistance management strategies must, however, be implemented before the resistance gene becomes common and stable in the population; otherwise, the resistant gene will not recede even if use of the insecticide causing selection pressure is discontinued. As demonstrated by a study of blowflies by McKenzie and Whitten in 1982 (3), fitness cost is not an intrinsic property of the gene. Therefore, if that gene is allowed sufficient time to become common in a population, the rest of the genome will adapt to incorporate it without a significant fitness cost. At this point, even if the selection pressure is removed, the resistance gene will remain in the population."

If this was the case it would probably favour rotations over sequences.

Do you have any thoughts ? Do you think it likely that costs will decline over time at high resistance frequencies ?

I looked at the 1982 paper (attached) and I don't think it supports the idea that costs decline over time (but it does seem to have data on relative fitness of RR,RS,SS that we could use in the WoS paper, gives evidence for dominance changing over time and could be used to derive our model parameters.)

Andy

Thanks Andy

That results should put the cat amongst the pigeons i.e. little difference between sequences and rotations.

There are examples of “compensatory” mutations offsetting the fitness cost of a mutation. Often from bacteria where populations can be very high and the population is asexual so that the mutation and modifier are always linked. It is less well studied in larger animals. There may be examples from IR but I am unaware. Often it is a post hoc explanation plucked out the air to explain why there is no apparent fitness effect i.e. “there were fitness effect, but they must have been modified”.

I like the 1982 paper… I have put it in the WoS folder with a note to discuss it.

25/1/2018

noticed that rotation plots don't have a log y axis where resistance plots do. There is a logy option in run_rot() should probably set that TRUE by default. yes fixed

Rowland(1981) has some good stuff on costs of resistance and selection coefficients and putting these into a model and testing against field data :

Withdrawal of dieldrin in two regions of India led to the frequency of SS increasing from 0.31 to 0.62 within 7 months (Rao et al., 1960) and from 0.11 to 0.85 within 9 months (Bhatia & Deobhankar, 1963). In Fig. 9 the computed reversions in An.garnbiae and Anstephemi have been redrawn, expressed in terms of SS rather than R frequency. Superimposed on these are the An. culicifacies field reversions. Assuming that one generation takes 1 month (cf. Curtis et al., 1978), the computer predictions seem to be consistent with the field reversions.

Having established why and to what extent resistance should revert when cyclodiene usage is ended, the question remains could yHCH be effectively re-introduced at a later date as part of a insecticide-rotation resistance management strategy? Because RR is the most deleterious genotype, the rate of reversion should be rapid initially. Only when R is predominantly carried by heterozygotes would the rate of reversion begin to slow down. For this reason resistance cannot be expected to revert completely, but it should nevertheless revert sufficiently to warrant the useful re-introduction of yHCH for several generations at a time. An unresolved question, potentially damagingto the long-term future of the rotation strategy,is would modifiers of fitness evolve if cyclodieneusage was prolonged?

It seems, then, that rotation strategies might prolong dieldrin efficacy indefinitely, not just against anophelines but against many species of insect.

[because the combination of selection coefficient, cost and dominance means that resistance doesn't rise too fast when used and declines when not in use].

What value of cost did Rowland1991 use ? Seems that he used high costs of 0.46 and 0.35, see table 8.

Rowland2013 data on hut trials of actellic - pirimphos-methyl - an OP. Includes monthly change and at diff concentrations so could be useful for WoS, but doesn't have resistant & susceptible (does have resistance frequencies).

Matowo2015 selection of insecticide resistance, seems to show changing insecticide resistance in response to interventions. but only shows mortality change not resistance frequency change.

Kim2016 interesting risk paper : The main policy implication of this analysis is that evolutionary speed of insecticide resistance and maximal larviciding efficacy strongly affect the preferred actions implied by our decision analysis tool and that these parameters jointly offer the largest potential benefits from reduced uncertainty

26/1/2018 new UI resistrot, similar name to resistmixseq

1/2/2018 I want to create a figure that shows effect of untreated area, refugia Can I do this in ggplot2, or do it in powerpoint ? see : resistance-dispersal-dilution-figure.pptx

How to create a greater effect of dispersal, reduce the coverage.

run_rot(n_insecticides = 1) Error in RAF1gen[, , "intervention"] : incorrect number of dimensions Called from: rot_migrate(RAF[, , , gen], migration = migration, coverage = coverage)

but this is OK run_rot(n_insecticides = 1, coverage=1)

I suspect problem is with dropping dimensions in rot_migrate when there's just one insecticide For now just force n_insecticides > 1 from the UI

~ c 15 generations per year for mosquitos, 25 day generation time. ~ 3 years = 45 generations ~ https://www.cdc.gov/malaria/about/biology/mosquitoes/

14/3/2018 re-red this Sternberg & Thomas 2017 Insights from agriculture for the management of insecticide resistance in disease vectors Because of the situation-specific nature of insecticide resistance management plans, public health may have an advantage over agriculture. Agricultural pest management is concerned with hundreds of species, and in many cases, detailed data are not available for each pest species (and potential nontarget organisms). Compared to agriculture, public health is concerned with a small number of insect species. This should enable a strongly data-driven approach to selecting resistance management strategies, but in reality, very little is currently known.

• The success of an IRM strategy is highly dependent on ecological and biological specifics; there is not yet a solid evidence base of recommending one strategy over another. • If control is achieved via the use of multiple tools with diverse modes of action, not only will selection for insecticide resistance be reduced but also, the impact of resistance will likely be less severe. The pending resistance crisis creates an urgent need to develop and implement integrated, multitactic IVM strategies that parallel IPM in agriculture.

To our knowledge, there is only one field trial aimed specifically at testing resistance management strategies against malaria mosquitoes (Penilla et al., 2007). Perhaps worryingly, this trial revealed no clear benefit of either insecticide rotations or mosaics in slowing the spread of resistance when compared to monotherapy. However, these results do not necessarily mean that rotations or mosaics cannot be effective, but rather highlights a critical gap in understanding when, how, and why such strategies might work.

15/3/2018 back from Peru and Geneva ~ get trial run of sensitivity analysis going ~ sensi_an_rotations1.Rmd

~ now have 500 runs working, max gens up to 400 ~ seems that generations are identical for sequence and rotations ~ results saved in df_res_all_rot.rda

Setting coverage to <1 (0.5) and migration to 0.01 does result in rotations lasting longer than sequences in some scenarios. Particularly shorter ones. Greatest slope is when rot_interval set to 5, one scenario gives 300 for sequence and 350 for rotation.

So maybe difference is created by the changeover, and more frequent changeovers create greater difference.

Also when coverage set to 1 again and cost > 0 there was 1 scenario with a very short rot-interval where the rotation lasts fair bit longer than the sequence.

scenario 79 expo_hi 0.79 eff 0.72 rot_interval 0

tot gens deployed under freq 0.5 = 199

scenario 79 expo_hi 0.79 eff 0.72 rot_interval 6

tot gens deployed under freq 0.5 = 251

run_rot(rot_interval=0, max_gen = 300, migration=0.01, expo_hi=0.8, eff=0.8, cost=0.1, coverage=1, add_gens_under50=TRUE) run_rot(rot_interval=6, max_gen = 300, migration=0.01, expo_hi=0.8, eff=0.8, cost=0.1, coverage=1, add_gens_under50=TRUE)

It could perhaps be just a rounding thing, because the number of generations is an integer there will always be some overshooting, when you have more switches there is the potential for more overshooting and the runs will last for longer.

~ checking on coverage=0, currently resistance frequencies increase when they shouldn't do run_rot(rot_interval=0, max_gen = 300, migration=0.01, expo_hi=0.8, eff=0.8, cost=0.1, coverage=0, add_gens_under50=TRUE)

in run_rot()

coverage is passed to : set_start_freqs() but only to determine whether the

RAF 3rd dimension site is set to : c('intervention','refugia') or just 'intervention'

actually if coverage=0 perhaps dimension should just be set to refuge.

So maybe it doesn't point to an error that resistance rises when coverage=0 just a property of the model structure.

16/3/18 ~ code condensing in run_rot()

~ added stop in run_rot() for if coverage=0

~ looking at comparing two scenarios (rot & seq) on same plot

I tried this before with : df_resanother in rot_plots. It puts on the resistance curves for the 2 scenarios, but doesn't show the insecticide use for each dfseq <- run_rot(n_insecticides=2, rot_interval=0, expo_hi=0.5, eff=0.5, cost=0.01, coverage=1) dfrot <- run_rot(n_insecticides=2, rot_interval=15, expo_hi=0.5, eff=0.5, cost=0.01, coverage=1) rot_plot_resistance(dfrot, df_resanother=dfseq, plot_refuge=FALSE)

dfseq <- run_rot(rot_interval=0, max_gen = 300, migration=0.01, expo_hi=0.8, eff=0.8, cost=0.1, coverage=1) dfrot <- run_rot(rot_interval=6, max_gen = 300, migration=0.01, expo_hi=0.8, eff=0.8, cost=0.1, coverage=1) rot_plot_resistance(dfrot, df_resanother=dfseq, plot_refuge=FALSE)

I have improved plotting of 2 scenarios together. Seems that overshooting problem occurs in rotations. Possibly because rotations are not stopped when the threshold is reached. See dfrot from prev e.g. whichj keeps going up to 0.62 at gen 127. This may not matter if the counting of gens only includes those < 0.5. The other option is to stop rotations early,m which I thought I was already doing with exit_rot = TRUE ??

in insecticide_check() exit_rot should prompt exit from rotation when 0.5 threshold is exceeded.

dfrot2 <- run_rot(rot_interval=6, max_gen = 300, migration=0.01, expo_hi=0.8, eff=0.8, cost=0.1, coverage=1, diagnostics=TRUE)

Aha! It could be a problem with min_rwr_interval, which is 5 by default and so may force staying within rotations and thus cause overshooting of threshold.

But i thought that only gens under thresh are counted ???

19/3/2018 on train to liverpool

trying to understand overshooting issue

but doesen't look like much time over the 0.5 thresh for the rotations can explain why they take longer. The short rotations only get to 0.5 towards the end.

Can I refine this example to understand what may be happening ?

turn off log ?

2 insecticides dfseq <- run_rot(rot_interval=0, max_gen = 70, migration=0.01, expo_hi=0.8, eff=0.8, cost=0.1, coverage=1, n_insecticides=2) dfrot <- run_rot(rot_interval=6, max_gen = 70, migration=0.01, expo_hi=0.8, eff=0.8, cost=0.1, coverage=1, n_insecticides=2) rot_plot_resistance(dfrot, df_resanother=dfseq, plot_refuge=FALSE, logy=FALSE) gens_under_thresh(dfrot) gens_under_thresh(dfseq)

This does go over the 0.5 threshold ? can I count this example in sensi_an_rotations ?

gens_under_thresh(dfrot) : 49 gens_under_thresh(dfseq) : 46

When there is no cost they are identical (39) : dfseq <- run_rot(rot_interval=0, max_gen = 70, migration=0.01, expo_hi=0.8, eff=0.8, cost=0, coverage=1, n_insecticides=2) dfrot <- run_rot(rot_interval=6, max_gen = 70, migration=0.01, expo_hi=0.8, eff=0.8, cost=0, coverage=1, n_insecticides=2) rot_plot_resistance(dfrot, df_resanother=dfseq, plot_refuge=FALSE, logy=FALSE) gens_under_thresh(dfrot) gens_under_thresh(dfseq)

but they still look diff in plots, i think because stoppping criteria are different in plots for seq & rot.

weirdly when I went to 3 insecticides, & cost 0.1

gen calc gives same result for each at 69 and plots finish at same point ?

aha! that was because max_gen=70 so max result could be 69 i.e. no insectidides in either scenario reached threshold

dfseq <- run_rot(rot_interval=0, max_gen = 100, migration=0.01, expo_hi=0.8, eff=0.8, cost=0.1, coverage=1, n_insecticides=3) dfrot <- run_rot(rot_interval=6, max_gen = 100, migration=0.01, expo_hi=0.8, eff=0.8, cost=0.1, coverage=1, n_insecticides=3) rot_plot_resistance(dfrot, df_resanother=dfseq, plot_refuge=FALSE, logy=FALSE) gens_under_thresh(dfrot) gens_under_thresh(dfseq)

1 20-20 = 0 2 insecticides 46 - 49 = 6 3 insecticides difference is 83-79 = 4 4 : 147 - 124 = 23 5 : 218 - 200 = 18 6 : 437 - 391 = 46

? no obvious pattern

done ~ write function to count gens under thresh done ~ add thresh line to plot done ~ add warning to run_rot() for if threshold not reached in time

exit_rot seems not quite to be doing what I expect. e.g. if I set to just 1 insecticide sequence keeps going longer than rot (although this is a somewhat weird example)

Aha! it is min_rwr_interval that did this. It stops a change early in a rotation but this can stop the simulation from ending.

(however it shouldn't effect the gens_under_threshold)

Try seeing if the gens_under_thresh are the same when min_rwr_interval is disabled.

(also note that min_rwr_interval isn't really sposed to apply to rotations, partly a result of the structure )

setting min_rwr_interval to 1 now results if diff of 86-77 = 9 for the 3 insecticides ??

20/3/18 Go back to the scenarios that was very different in sensi_an_rotations_test num79

To enable side by side comparison

library(gridExtra) ggrot <- rot_plot_resistance(dfrot, plot_refuge=FALSE, logy=FALSE) ggseq <- rot_plot_resistance(dfseq, plot_refuge=FALSE, logy=FALSE) grid.arrange(ggrot, ggseq, ncol = 1)

I think maybe need to go back to two insecticides to try to find small difference :

ni <- 2 cost <- 0.1 gens <- 70 dfseq <- run_rot(rot_interval=0, max_gen=gens, migration=0.01, expo_hi=0.8, eff=0.8, cost=cost, coverage=1, n_insecticides=ni) dfrot <- run_rot(rot_interval=6, max_gen=gens, migration=0.01, expo_hi=0.8, eff=0.8, cost=cost, coverage=1, n_insecticides=ni) ggrot <- rot_plot_resistance(dfrot, plot_refuge=FALSE, logy=FALSE) ggseq <- rot_plot_resistance(dfseq, plot_refuge=FALSE, logy=FALSE) grid.arrange(ggrot, ggseq, ncol = 1) gens_under_thresh(dfrot) gens_under_thresh(dfseq)

rot:48, seq:46

Note that in the above plot other strange issues, which may be plotting probs. ~ insecticide use boxes don't join up, seem to be gaps ~ resistance keeps increasing after insecticide box has finished ~ first rotation seems shorter than later, because insecticide is classed as NA in firs tstep, is that right ?

Difference is for i2 which has 4 rots of 6, 24, only reaches threshold at end of final rot In seq, 12 last for 2 generations less.

Trying too look at plot in more detail dfrot2 <- run_rot(rot_interval=6, max_gen=20, migration=0.01, expo_hi=0.8, eff=0.8, cost=cost, coverage=1, n_insecticides=ni)

Not sure if I'm going to be able to fix the ribbons joining ? can I +1 to the x upper limit ?

to_plot_in_use <- 1

dfrot2 <- dfrot2 %>% mutate( i_in_use = ifelse(resist_gene==paste0('insecticide',insecticide), to_plot_in_use, NA))

gg <- ggplot( dfrot2, aes_string(x='generation',y='resistance')) + geom_line( alpha=0.5, lwd=1, colour='red3' ) + facet_wrap('resist_gene', ncol=1) + geom_ribbon( aes_string(x='generation', ymin=0, ymax='i_in_use'))

plot(gg)

I think ribbons not meeting and curve increasing after are just plotting issues that I don't know how to fix.

If a rotation is 10 generations, the first finishes at 10 and the next is shown as starting at 11. The effect of the first continues until 11. Setting i_in_use to 1 for the following generation would fix the plotting issue.

Do I need to try to add a generation 0 to store starting conditions ? NO

Generations The conditions within a generation influence resistance which is stored in the following generation. Therefore it's good to set insecticide for time 1, and this is consistent with later generations.

So I've fixed saving of insecticide for t1

Still have difference between 24 & 22 for i2 in above example

add in_use columns

dfrot <- dfrot %>% mutate( i_in_use = ifelse(resist_gene==paste0('insecticide',insecticide), 1, NA)) dfseq <- dfseq %>% mutate( i_in_use = ifelse(resist_gene==paste0('insecticide',insecticide), 1, NA))

remove NA refuge rows

dfrot <- dfrot[dfrot$active_or_refuge=='active',] dfseq <- dfseq[dfseq$active_or_refuge=='active',]

write_csv(dfrot,'rot.csv') write_csv(dfseq,'seq.csv')

see : rot_seq_compare_march2018.xlsx

With cost it seems that rotation may be able to last longer than a sequence.

Could this be because decline in resistance freq is more rapid at higher frequencies.

Both increases and decreases in resistance frequency are higher at higher resistance frequencies.

But does this explain ?

Is it something to do with more frequent changes ?

Could it be because : ~ frequency changes are bigger as frequencies get higher, ~ allowing in a sequence to reach high frequencies straight away makes it slightly quicker than see-sawing at lower frequencies ~ is their some maths to prove this ? ~ ask Ian ~ can I print an example ?

Ian said that at higher dominance the difference in change between low & high frequencies will be greater, so might expect the difference between rot and seq to be greater. Try it and see.

linputs <- list()

linputs$n_insecticides <- 2 linputs$coverage <- 1

linputs$migration <- 0.1

linputs$cost <- 0.1 linputs$max_gen <- 200 linputs$expo_hi <- 0.6 linputs$eff <- 0.6 linputs$dom_sel <- 1 linputs$dom_cos <- 1

linputs$rot_interval=0 dfseq <- do.call(run_rot, linputs) linputs$rot_interval=6 dfrot <- do.call(run_rot, linputs)

ggrot <- rot_plot_resistance(dfrot, plot_refuge=FALSE, logy=FALSE) ggseq <- rot_plot_resistance(dfseq, plot_refuge=FALSE, logy=FALSE) library(gridExtra) grid.arrange(ggrot, ggseq, ncol=1) gens_under_thresh(dfrot) gens_under_thresh(dfseq)

with higher dominance there is only 1 generation difference (31:30) may be partly because it is faster

when I reduced exposure and effectiveness then a largish difference between rot & seq did return.

stackoverflow says this is a cygwin memory error and rebooting should fix it. eek 0 [main] us 0 init_cheap: VirtualAlloc pointer is null, Win32 error 487 AllocationBase 0x0, BaseAddress 0x68570000, RegionSize 0xE0000, State 0x10000 C:\Program Files (x86)\Git\bin\ssh.exe: *** Couldn't reserve space for cygwin's heap, Win32 error 0 fatal: Could not read from remote repository.

Just rebooted and that fixed it.

Seems that rotations favoured by higher costs, higher dominance of cost, lower start freq (longer to act), shorter rot interval (more changes)

3/4/18 ~ set up initial sensi analysis of 1000 runs ~ 200 runs take < 30 mins

~ seems that a lot of runs don't reach resistance thresholds by 400 generations, and this is with start-freqs set to 0.01, if I set lower start freqs then even more runs will be excluded.

done ~ add ranges for variable params done ~ copy across grid.arrange from phaps-rotations-better

around 79% of runs failing to make thresholds

upped max gens to 500 and the lower exposure val to try to improve number making thresholds

4/4/2018 ooo with new ranges getting a few runs where there is a big advantage to rotation. Seems to be assoc with high dom_cos > 0.8 Only c 1% of scenarios, 2 out of 200

scenario 137 sum generations deployed insecticide under thresh rotation: 500 sequence: 258

run_rot(n_insecticides=3, cost=0.026, start_freqs=0.032, rot_interval=10, eff=0.5, dom_sel=0.36, dom_cos=0.82, rr=0.32, expo_hi=0.73, coverage=1, migration=0, max_gen=500, min_rwr_interval=5)

scenario 143 sum generations deployed insecticide under thresh rotation: 500 sequence: 126

run_rot(n_insecticides=2, cost=0.049, start_freqs=0.064, rot_interval=0, eff=0.71, dom_sel=0.089, dom_cos=0.87, rr=0.9, expo_hi=0.65, coverage=1, migration=0, max_gen=500, min_rwr_interval=5)

In 137 even the rounding makes the rotation last slightly less long (e.g. c400 rather than 500)

start writing code to calc difference between sequence and rotation ready to put that into PRCC.

currently i have repeated following rows for sequence and rotation with identical input params

i probably want to replace this with one row and columns for gens_rot and gens_seq

new version one row per scenario :

advantage to rotation despite dom_cos=0.23, cost=0.02

scenario 153 sum generations deployed insecticide under thresh rotation: 406 sequence: 289

run_rot(n_insecticides=4, cost=0.023, start_freqs=0.019, rot_interval=28, eff=0.92, dom_sel=0.019, dom_cos=0.23, rr=0.63, expo_hi=0.75, coverage=1, migration=0, max_gen=500, min_rwr_interval=5)

Looking at some more of literature :

From Mallet 1989 : Since a naive strategy - of using insecticides in turn until resistance evolves to each - gives almost the same total control period as a planned rotation23, rotations can hardly even be defined as management. However, in some cases rotations may be the only practical way of limiting resistance to a single favoured compound 23 27.

23 is curtis 1987 a book chapter I couldn't find.

Curtis 1993 has some good quotes and references about rotations. It covers the oncheriasis control program which is taken by some to be evidence in favour of rotations Fig 3 & 4 is lab data comparing sequence & rotation, very similar to our plots, I should definitely cite it.

"A programme such as the OCP differs from these simple models in that more than two insecticides are available for use and a given compound can therefore be out of use for more generations than it is in use."

I hadn't thought of that. Having more than 2 insecticides allows each to be out of use for more time than it is in use. So if the increase in resistance while in use was the same as the decrease when not in use then a rotation could be sustainable in the long term. However it is likely that resistance will increase faster in presence than it will decrease in absence (I think Curtis says that too).

Curtis cites how his 1987 work has already included effect of modifier genes, and shows it here in Fig 6. "Supporters of prearranged rotations usually tacitly or explicitly assume that the general fitness due to resistance genes will improve if they are exposed to selection for many generations. Curtis (1987) simulated the behaviour of a modifer gene which could have this effect and found that it would only have an appreciable influence if one allowed resistance gene frequencies to become very high before switching insecticides.""

"One possible objection to the reactive system is that it requires an efficient and comprehensive system for testing for resistance to avoid ‘overshooting’ and getting stuck as in the dotted line in Fig. 6. However, running a prearranged rotation without a monitoring system would mean that one would have no way of knowing if the system was still working."

on p16 there is also interesting consideration of what expected in a mixture if one insecticide decays before another.

It's negative about rotations : 'Except where there are seasonal/environmental reasons for a pre-arranged rotation (as in the OCP), we do not find that such rotations are likely to have any advantage over well-managed reactive rotation, and we do not consider that the latter can validly be called a 'strategy' for resistance management. The use of a suitable insecticide mixture seems more promising in the light of the limited data so far available on resistance in relation to insecticide impregnated bednets.'

I would like to get the Curtis 1987 paper Curtis CF. 1987. Genetic aspects of selection for resistance. In Ford MG, Holloman DW, Khambay BPS, Sawicki RM, eds. Combating resistance to xenobiotics. Chichester, U.K.: Ellis Horwood, 15e161.

but couldn't find and probably adds little

5/4/2018 added progress counter to sensi_an_rotations after 2000 runs failed to finish overnight also now save the output every 200 runs

100 scenarios take < 10 mins 2000 scenarios should take 20*10 mins, 200 mins < 2 hours but seems that RAM fills up that may slow it down

got 1400 runs, PC ground to a halt.

No very clear pattern for high difference runs. Small positice effect with cost & dom_cos, small -ve effect of dom_sel

9/4/2018 turned off plotting in sensi_an_rotations 10,000 runs took c 100 minutes (1hr20)

27/4/2018 First draft of paper3

trying to get references to work without having to copy the .bib file from resistance

C:/rsprojects/resistance/inst/documents/paper2/

Getting filter pandoc-citeproc' had status 83 Execution halted https://stackoverflow.com/questions/36826916/pandoc-citeproc-error-83-with-rmarkdown-file

removed all '-' from citation keys in mendeley, didn't fix

Temporary solution of cutting out all refs except 2 from the .bib file fixed it. Need to go throught the bib file and look for error. Can I reduce size of the bib file to just the IR papers ? What did I do for ecocrop ?

30/4/18

100 runs with no cost all had identical times for rotations and sequences

from Rmd got error when tried 1000 worked from console

With no cost and no migration, majority of runs the same. Just 3/1000 seq better than rot by < 3 generations. tst <- filter(df_res_all, rot_minus_seq!=0)

Now trying 1000 runs with coverage and dispersal set between 0.1 and 0.9

Seems that less benefit to rotations with dispersal rather than costs.

2/5/18 exchanging emails with Ian

Thanks Ian, I'm getting closer to understanding. I'm just trying to understand why the lines cross, which suggests to me that resistance wouldn't be selected for at low frequencies ?

This is a simple text description of how I'm starting to understand. Is it true ? (we can work on improving the wording later).

Resistance increases when an insecticide is in use and, if there are resistance costs, decreases when it is not.

If the difference between this increase and decrease over periods of use and non-use is positive then resistance will keep increasing.

If the positive difference between the increase and decrease is lower at lower frequencies then there is an advantage to short rotation intervals that keep frequencies low.

In this latter scenario sequential use leads to higher frequencies where increases during periods of use are more likely to outweigh decreases during periods of non-use.

done ~ Can I change histograms of count of rot advantage to have different bins, we're only interested in whether difference is 0, moderate or substantial.

Trying to recreate scenario 2167 after removal of non threshold runs :

scenario <- 2167 linputs <- as.list( df_cost[scenario, 1:15] ) dfres <- do.call(run_rot, linputs) ggrot <- rot_plot_resistance(dfres, plot_refuge=FALSE, logy=TRUE, plot=FALSE) linputs$rot_interval <- 0 dfres <- do.call(run_rot, linputs) ggseq <- rot_plot_resistance(dfres, plot_refuge=FALSE, logy=TRUE, plot=FALSE) grid.arrange(ggrot, ggseq, ncol=1)

linputs$dom_cos <- 0.8 x11()

linputs$rot_interval <- 10 dfres <- do.call(run_rot, linputs) ggrot <- rot_plot_resistance(dfres, plot_refuge=FALSE, logy=TRUE, plot=FALSE) linputs$rot_interval <- 0 dfres <- do.call(run_rot, linputs) ggseq <- rot_plot_resistance(dfres, plot_refuge=FALSE, logy=TRUE, plot=FALSE) grid.arrange(ggrot, ggseq, ncol=1)

setting params for 2167

linputs <- list() linputs$n_insecticides <- 3 linputs$coverage <- 1 linputs$migration <- 0 linputs$cost <- 0.0582 linputs$max_gen <- 500 linputs$expo_hi <- 0.64 linputs$eff <- 0.96 linputs$rr <- 0.78 linputs$dom_sel <- 0.06 linputs$dom_cos <- 0.7 linputs$start_freqs <- 0.0731 linputs$male_expo_prop <- 0.167 linputs$plot <- FALSE

this shows a small tweek of the params leads to the sequence lasting forever too

and I don't know which input change causes the seq to last forever ?

Aha seems that I hadn't included rr=0.78, rr default is 0.5,

with lower rr less selection for r and sequence does better

linputs <- list() linputs$n_insecticides <- 3 linputs$coverage <- 1 linputs$migration <- 0 linputs$cost <- 0.058 linputs$max_gen <- 500 linputs$expo_hi <- 0.64 linputs$eff <- 0.96 linputs$rr <- 0.5 #0.78 linputs$dom_sel <- 0.06 linputs$dom_cos <- 0.7 linputs$male_expo_prop <- 0.167 linputs$start_freqs <- 0.0731 linputs$plot <- FALSE

3/5/18 Targets for today a) check how stable the rotations advantage is b) understand and display the mechanism

But I went off that track because I got more to thinking we should be trying to identify circumstances where both strategies 'succeed'. However I suspect that 'success' here is due to low effectiveness and exposure which lead to slow evolution but possible cause low control too.

Looking at trees to differentiate between success, failure and rot_only. Not very helpful.

Looking at input effects on n generations, they go the way I would expect, but may have to be part of the story.

The violin plots for each input var for rot_seq, rot, neither look pretty good.

e.g. for cost they show that successful runs tend to have higher costs, those failing for both have lower costs. Those that work just for rotations tend to have lower costs.

8/5/18 Getting violin plots working.

Maybe a more useful question, rather than when are rotations or sequence better is when are both rotations and sequences sustainable. i.e. under which combination of conditions do both strategies last 'forever' ?

In that case the runs that I have been excluding become the useful ones.

Because actually, finding the runs that are different for rotations and sequences after excluding those that 'succeed' for both is actually trying to find those that 'fail' for sequences.

I think there could be a boundary area in between those that 'succeed' for both and those that 'fail' for both, where some 'succeed' only for rotations.

9/5/18 Leaving in all runs now

Something weird has happened to fig1 in the paper3 ?? because of a > in fig legend

Also because we are leaving all runs in, the indices of the runs with large rot advantage are now wrong ! I should have recoded a unique id for runs.

ARg, now all plots have stopped appearing in doc and I don't know why ????? It was a < symbol in the MS text, could tell because the text stopped at that point.

10/5/18 Testing setting width of plots for publication at 300 dpi into word doc.

Solution : just use fig.height & width (out.height etc not supported for word) This makes the tiffs correct. Word foten resizes them. To get right in word, right click set size & position, just click reset to go back to 100%.

improved histogram figures

got Ians fig Y working in freq_fitness_simple() & added facetted grid plot to the paper.

11/5/18 - improved violin plots added to pap3 - improving legends

Finally grasping the nettle of insecticide switching.

I want to stop switching to an insecticide that was used < X generations ago.

To do that store an array which has a counter for each insecticide which records how long since it was last used. Then just add to that counter in each generation for each insecticide not in use. Don't allow an insecticide to be used that is below a user-defined level.

in insecticide_switch() BE CAREFUL!

min_gens_switch_back minimum num gens before can switch back to an insecticide

I set the default in run_rot to 10 generations. (I can leave min_rwr_interval in for now, when it is set to 1 it has no effect).

Note that min_rwr_interval was assessed in insecticide_check() i.e. it stopped an insecticide from being stopped within a certain interval

Which actually maybe wasn't so bad ?? Maybe I could stick with it ?? i.e. that simulated not being able to withdraw an insecticide which has already been applied (which sounds pretty good now I think about it ...)

Either of these could be OK. I had it set to 5 generations before which was half a year.

done ~ add to methods about min rwr interval done ~ add about keeping resistance allele freq at starting level done ~ add to fig Y legend Ians caveats

14/5/2018 ~ Paper3, decsribing results and adding to methods clarifiying model constraints.

Before sending to Ian. ~ resize non square figs, make sure one per page

looking at how I made smaller reference bib file for uea work.

Seems I can export a single non-synced bibtex file

https://www.imperial.ac.uk/media/imperial-college/administration-and-support-services/library/public/LaTeX-and-Mendeley-sep-2017.pdf 1. Exporting a BibTeX (.bib) file from your Mendeley library • Select the records for export • From top menus, go to File>Export… • In Popup box, choose a file name (no spaces); ensure file type is set as BibTeX • Save the .bib file in the same folder as the LaTeX document

df_res_all_rot_nocost_nomig_10000_mgs5.rda

mgs5 to indicate min_gens_switch_back=5

Aha! all these big sequence advantages (well those >50%) were due to artfefact when rotation interval=5 and num insecticides=2 rotation couldn't return to first insecticide because it was only 5 generations ago !!

17/5/2018

New papers to cite in rotations paper :

Hackett (2016) https://www.ncbi.nlm.nih.gov/pmc/articles/PMC5026168/ interesting model of resistance to transgenic crops with similar results on the importance of dominance of resistance and cost : "the most significant constraint upon model behaviour is the dominance of the resistant phenotype, with greater dominance decreasing Bt usage. For a given level of dominance of resistance, the type and dominance of any fitness costs may relax these constraints, allowing for larger areas of Bt crop [smaller refuges to control resistance evolution]"

The R code is here and looks quite accessible : https://datadryad.org//resource/doi:10.5061/dryad.4s405

Allen(2017) Reversing resistance: different routes and common themes across pathogens http://rspb.royalsocietypublishing.org/content/284/1863/20171619

Mostly applied to microbes & fungi, some good general points, probably not useful to us. Does talk about different types of reversion and the importance of cost.

My sentence from an example about negative cross resistance : Negative cross resistance meant that two fungicides were used successfully on Botrytis cinerea in European vineyards for 7 years ending in 1994 when a new mutation appeared that gave resistance to both. [87]

Russell PE. 1995 Fungicide resistance—occurrence and management. J. Agric. Sci. 124, 317–323. (doi:10.1017/S0021859600073275)

So negative cross resistance can work for a while but is susceptible to new mutations giving resistance to both ingredients.

Alphey & Bonsal 2017 good accessible description of high dose-refuge strategy : Genetics‐based methods for agricultural insect pest management (but the paper itself is about sterile insect technique)

For example, modelling the effects of larval competition and exploring late‐acting lethal phenotypes in mosquitoes predicted that this could be substantially more effective at population control than an early‐acting (e.g. embryonic) lethality or radiation‐induced sterility (Atkinson et al., 2007; Phuc et al., 2007; Alphey & Bonsall, 2014a). Indeed, if density‐dependent juvenile competition were over‐compensatory, genetic lethality that occurred at an earlier stage, thereby freeing survivors from regulation by intense competition, could push adult insect numbers higher than in the natural uncontrolled population.

In terms of the genetics, consider a resistant allele r, which is initially rare. The dominant allele S is susceptible to Bt. If the high dose assumption is achieved, only some rr individuals can survive a full life cycle on transgenic plants and emerge as adults to mate. In the refuge, most emerging adults will be susceptible SS, especially if the r allele has fitness costs in the absence of the Bt toxin. If the refuge is located so that the two subpopulations are well‐mixed, most resistant rr survivors will mate with susceptible SS insects from the refuge. Their resulting Sr progeny cannot survive on the Bt crops and so will not pass on the resistant allele to future generations. Unfortunately, a few insect species, such as the economically important pests Helicoverpa zea (Boddie) and Helicoverpa armigera (Hübner) (both species, confusingly, known as both bollworm and corn earworm), have been identified in the past as ‘moderate dose pests’, where the toxins were unable to kill all heterozygotes (EPA, 1998; Tabashnik et al., 2008). Even where its main assumptions appear to hold, the high‐dose/refuge strategy is predicted only to delay resistance and, after two decades of commercially grown Bt cotton and Bt corn (maize), some field‐evolved resistance has now been observed and reported in a variety of insect species. In some cases, this has already led to reduced efficacy of crops, or even crop failures. A comprehensive review is provided by Tabashnik et al. (2013).

In economic terms, the high‐dose/refuge strategy is an inter‐temporal trade‐off, sacrificing current value to retain value generated in future (Frisvold & Reeves, 2008). Refuge plants will be damaged, which reduces their yield if they are crops, or reduces the area allotted to crop production if alternative host plants are used. This damage or lost production is tolerated, in return for prolonging the efficacy of the protection afforded by the Bt crops. In principle, a conventional crop refuge could be sprayed with another pesticide (one with no cross‐resistance between its active ingredient and Bt), although doing so reduces its effectiveness as a refuge.

Aha and a justification for using resistance threshold of 0.5 ! Generally, a frequency of 0.5 is an indicator of a serious problem because resistance frequency increases rapidly once r becomes the more common allele (Carrière & Tabashnik, 2001).

Carrière & Tabashnik, 2001. Modelling refuges in Bt corn Good paper to cite in rotations paper. They found some similar things to us, e.g. resistance not increasing when dominance of resistance low and dominance of cost high.

Carriere 2018 good one to include, recent data on fitness costs showing that there are highish fitness costs but that they are almost completely recessive. Genotype-specific fitness cost of resistance to Bt

18/5/2018

Glunt 2018 Impact of temperature on insecticide toxicity. Relevant to window of selection. In arabiensis & funestus with deltamethrin, generally mortalities higher when temperatures either higher or lower than standard. The difference between resistant and susceptible mortality is generally less, suggesting that the window of selction may be lower outside of standard lab conditions. Bendiocarb mortality increased with temperature and unsurprisingly was similar between strains susceptible and resistant to pyrethroids.

21/5/2018 Gould et al. Science paper. Wicked evolution: Can we address the sociobiological dilemma of pesticide resistance?

If we are to address this recalcitrant issue of pesticide resistance, we must treat it as a “wicked problem,” in the sense that there are social, economic, and biological uncertainties and complexities interacting in ways that decrease incentives for actions aimed at mitigation. Here, we summarize the interacting factors and conclude with a call for government support of ambitious landscape-level experiments to assess which pesticide use strategies decrease resistance risks.

In sexually reproducing weeds and insects, the rate of resistance evolution is strongly influenced by the relative fitness (dominant to recessive) of heterozygotes, and this sometimes depends on the dose of pesticide encountered in the field (Fig. 3). Thus, it is difficult (and controversial) to determine whether resistance is expected to evolve more rapidly to higher or lower application concentrations of a pesticide [e.g., 16, 17]. Even more complexity arises in attempts to predict resistance evolution when combinations of pesticides are applied (18, 19). Although the idea that such combinations will slow resistance evolution is theoretically controversial and lacks empirical support, mixtures are often recommended at the field level (15).

Model-based analysis has shown that if disease vector resistance to pyrethroids becomes widespread, cases of malaria averted with ITNs could decline by 40% (25). Coupled with the estimate that bednets averted more than 65 million clinical malaria cases in sub- Saharan Africa in 2015 (7), and assuming that this figure provides a lower bound for potential cases averted in subsequent years, this would imply around 26 million additional clinical cases of malaria per year as a result of widespread vector resistance. Assuming an approximate lower bound cost of illness of at least $10 per malaria episode (26), insecticide resistance could conservatively cost sub-Saharan Africa at least $260 million per year.

24/5/18 replotted PRCC for all scenarios

resistance cost comes out as a negative impact on benefit of rotations. This could be because higher cost runs tend to be equal because they both sequences and rotations don't reach resistance thresholds. Maybe the PRCC is confusing ?

25/5/18 While out for a run I thought that the constraint on allele frequencies not going below their starting values could have implications for difference between rotations & sequences. I did 10k cost simulations with the constraint removed to check. Either way I should include something about this in the discussion.

2/7/18 back to rotations after while on wos

try to knit rotations paper get : error in do.call ... object 'run_rot' not found Build & reload fixed it.

Also the citations aren't being picked up now ? Because I messed up export of library.bib from Mendeley, exporting just a single selected reference.

4/7/2018 looking at Hackett2016, downloaded R code which is <100 lines and the script works to produce the first figure.

I was thinking that this could be a simpler model that we could also use with the same inputs that we use currently. This could possibly allow us to look at popn size, and to have a quick assessment of when resistance will increase and decrease.

The growth of the pest population is described by a density-dependent selection model (Roughgarden 1971) which uses the two parameter density dependence function developed by Maynard Smith & Slatkin (1973).

But maybe the way it works isn't appropriate for mosquitoes :

Pests divided between refuge and toxic patch
Pests within the Bt patch are exposed to toxins and mortality dependent on genotype.
Pest mortality via intraspecific competition (incl. costs of resistance)
Surviving larvae mature and mate at random to generate the next generation.

In mosquitoes competition likely to occur at larval stage, maybe doesn't matter ?

9/7/18 mutation-seln balance

not quite working in run_rot() yet

run_rot(no_r_below_start=FALSE,no_r_below_mut = TRUE)

Ian said that mutation rate of 10e-9 is used in malaria studies. John Maynard Smith, Evolutionary Genetics p62 suggests mutation rate per generation per locus from drosophila of 10e-6.

Ian suggested that he expected these frequencies to be so low that using them as a lower limit will not effect simulations.

Checking ability to do runs where insecticides differ. Code at end of sensi_an_rotations1

Currently getting these warnings, & seemingly more rows than scenarios in df_res_all ... 1: In if (dom_cos == 0) { ... : the condition has length > 1 and only the first element will be used

fixed vectorising mutn_seln_bal.r

Running for insecticides with different inputs now works : dfr <- run_rot(n_insecticides=3, eff=c(1,0.6,0.4), rot_interval=0)

but I don't get one row per scenario in df_res_all from sensi_an_rotations

for 10 runs I get 31 rows because I'm getting 1 row per insecticide in each run where does that happen ? seems it happens in this line :

df_in_out <- as.data.frame(c(linputs, gens_rot=gens_rot, gens_seq=gens_seq))

when linputs contains some elements that are vectors (repeated for each insecticide) the resulting dataframe has one row per insecticide.

I would rather have one row per scenario and have all the vector elements saved to one record, perhaps as a comma delimited string ?

seems that as.data.frame is the culprit. what about tidyverse option ? no didn't work

this is the issue I have : as.data.frame(c(double=c(9,6), single=7))

converts to this : double1 9 double2 6 single 7

when I want something like double1 9,6 single 7

13/7/18 ltst <- list(d=c(9,6),s=7)

This nearly does what I want but loses names of columns paste0(ltst) [1] "c(9, 6)" "7"

cool this does it : lapply(ltst,toString) $d [1] "9, 6" $s [1] "7"

linputs2 <- lapply(linputs, toString) df_in_out <- as.data.frame(c(linputs2, gens_rot=gens_rot, gens_seq=gens_seq))

Seems that with variable insecticides there are more runs where sequences are better (i.e. highly negative rot_minus_seq) ? See : df_res_all_rot1000-diffI-lowstart.rda

17/7/18 To be able to plot specific scenarios with different insecticides I need to get the strings that I saved above back into vectors.

strsplit seems it should do the job :

tst <- as.character(linputs$rr) unlist(strsplit(tst,split=", ")) #note the space after the comma

30/7/18 continuing getting it to cope with multiple insecticides in sensi_an_rotaions

I'm trying to convert string representations of mutliple parameters (that were needed to save in an output file) back to vector versions so I can rerun function with them.

linputs2 <- lapply(linputs, function(x) ifelse(grepl(",",x), as.numeric(unlist(strsplit(as.character(x),split=", "))),x))

Currently I just get the first value

linputs$cost [1] "0.0209911216050386, 0.00808626485522836, 0.0159290118375793" linputs2$cost [1] 0.02099112

There was something a bit weird about the ifelse that meant just the first dimension was returned. This now works :

linputs <- lapply(linputs, function(x) if (!grepl(",",x)) x else as.numeric(unlist(strsplit(as.character(x),split=", "))))

Examples with different insecticides where sequences last much longer than rotations seem to be caused by one of the 3 insecticides lasting much longer in the sequence, in the rotation it gets to a stage where it is the only one remaining, and our changing algorithm doesn't allow switching back to it from itself.

~ added id to outputs in sensi_an_rotations to make things easier in future, e.g. to identify runs with particular behaviour. ~ also made more efficient by declaring size of results df at start of loop

The scenarios with the biggest differences seem to be where insecticide1 is the better one. Just want to check it isn't something else screwy ???

31/7/18 Getting this error pushing to github : 0 [main] us 0 init_cheap: VirtualAlloc pointer is null, Win32 error 487 AllocationBase 0x0, BaseAddress 0x68570000, RegionSize 0x340000, State 0x10000 C:\Program Files (x86)\Git\bin\ssh.exe: *** Couldn't reserve space for cygwin's heap, Win32 error 0 fatal: Could not read from remote repository.

Please make sure you have the correct access rightsand the repository exists.

updating RStudio didn't fix, had once before and rebooting fixed it then ...

rebooting didn't fix, this post suggested updating git to v2 from https://gitforwindows.org/ https://stackoverflow.com/a/31970708/1718356

git --version git version 1.9.5.msysgit.0

now I have 2.18.0.windows1, the installer sorted uninstalling the old version, and I didn't have to change anything in RStudio.

All rotation runs before were based on insecticides with equal effectiveness, resistance restoration, cost etc. These new runs allow differences between insecticides.

One thing that they showed was that a sequence could do much better than a rotation due to the logic of when to change an insecticide.

because just one insecticide remained below the resistance threshold and our logic didn't allow it to keep being used. These were generated when I allowed insecticides to differ. Below is the result for the same inputs of adding a condition to check all the other insecticides, and to allow us to keep using one if it is the last remaining. I have set the minimum rotation interval to 5 otherwise it keeps switching back to the other insecticides for a single generation. This gives more similar results for the rotation and sequence.

1/8/2018 Reading Curtis(1987) Genetic aspects of selection for resistance. Good short article. Summary : ~ rotations and sequences very similar ~ fitness modifiers unlikely to give benefit to rotations (unless resistance thresholds and costs very high) ~ there maybe operational reasons for favouring sequences

I should try to make sure our rotations paper is as easy to read as this.

It is a good justification for not using modifiers so I don't think we need to implement modifiers but if we did, this is how he describes what he did - although all the details not there. "assumed a single modifier of each major resistance gene. The fitness disdavantage of each major gene was eliminated when the modifier was homozygous and halved when it was heterozygous. The modifier itself caused slight fitness disadvantage and was initially maintained in equilibrium by a blance of origin by mutation and elimination by selection."

Other papers in the book: Denholm, Laboratory simulation of selection for resistance. Pyrethroid resistance in houseflies. Compares model results and lab experiments. Concludes that more field work needed before models can be operationally useful.

Dennehy, Decision making for managing pest resistance to pesticides. Brief possibly useful stuff on decision making criteria in agriculture.

Sawicki, chap9. Definition, detection and documentation of insecticide resistance.

Tried a solution to the diff insecticide problem. If no other insecticides are below resistance thresh and this one is then can keep using it in a rotation ...

good refs for rotations paper

resistance : malERA2017_updated-research-agenda-resistance-in-malaria.pdf

"What are the benefits of insecticide rotations, mixtures, or spatial mosaics of different compounds? What is the impact of adding nonpyrethroid IRS where LLINs are already deployed at high coverage and quality? When should new insecticides be adopted? What is the ideal rotation period or mosaic configuration? How many insecticide classes are needed for effective rotation or mosaic strategies? Despite the absence of data to answer these questions, some countries have already developed operational frameworks for resistance management that could be adopted by other programmes.""

The malERA Refresh Consultative Panel on Insecticide and Drug Resistance (2017) malERA: An updated research agenda for insecticide and drug resistance in malaria elimination and eradication. PLoS Med 14(11): e1002450. https://doi.org/10.1371/journal.pmed.1002450

Slater2016. Use of an individual-based simulation model to explore and evaluate potential insecticide resistance management strategies

12/11/2018 Coming back to this after a while before presenting to IVCC.

from UI : Warning: Error in raf_get: object 'RAF' not found but works for 2nd plot on right ?

must be from run_rot() I don't understand why sceanrio A fails and B doesn't ?

178: raf_get [C:/rsprojects/rotations/R/insecticide_check.r#70] 177: insecticide_check [C:/rsprojects/rotations/R/insecticide_check.r#70] 176: run_rot [C:/rsprojects/rotations/R/run_rot.r#325]

Also fails from the command line :

run_rot( n_insecticides = 3 , max_gen = 150 , start_freqs = 0.01 , rot_interval = 10 , eff = 0.8 , dom_sel = 0.5 , dom_cos = 0.5 , rr = 0.5 , expo_hi = 0.5 , coverage = 0.8 , migration = 0.01 , cost = 0.01 , logy = TRUE , no_r_below_start = TRUE , min_rwr_interval = 5 )

Error in raf_get(RAF, insecticide = -current_insecticide, sex = "f", gen = 1, : object 'RAF' not found

BUT now seems to be working from the command line ??? and is working from shiny ... weird.

this error then popped up again from shiny later on.

13/11/2018

preparing slides for IVCC talk next week

~ add a slide about decision rules : because just one insecticide remained below the resistance threshold and our logic didn't allow it to keep being used. Below is the result for the same inputs of adding a condition to check all the other insecticides, and to allow us to keep using one if it is the last remaining. I have set the minimum rotation interval to 5 otherwise it keeps switching back to the other insecticides for a single generation. This gives more similar results for the rotation and sequence.

23/11/18

to output script as pdf for Ian to print

library(knitr) stitch(script="R/run_rot.r")

stitch(script="R/run_rot.r", system.file("misc", "knitr-template.Rhtml", package="knitr")) stitch(script="R/rot_migrate.r", system.file("misc", "knitr-template.Rhtml", package="knitr"))

Then open html in a browser, print to pdf, set landscape and scale=150

5/12/18

For rotations sent Jim a scaled down version of the rotations sensitivity analysis code. That just runs the sensitivity analysis and Jim will have a go at parallelising it.

6/12/18

Check why git changes are thinking that whole file has changed. How did I fix that for rnaturalearth ??

GitHub suggests that you should make sure to only use \n as a newline character in git-handled repos. There's an option to auto-convert:

git config --global core.autocrlf true #global didn't work git config core.autocrlf true #just this project

Doing it for rotations project did fix the problem straight away. I didn't even have to restart RStudio.

6/2/19

before meeting with Jim to look at rotations code previous error in raf_get() seems OK now

13/2/2019

looking at saving tested mortality levels for the different insecticides over time (as a prelude to using it as a criteria for switching insecticides)

resistance frequency results are saved in df_results think I can create df_mort with the same structure and then put it through the same plotting function.

but df_results summaries are only calculated at the end so df_mortali couldn't allow mortality to be used as a switching criteria.

BUT also I need an equivalent of RAF[] (without the sex dimension, or maybe easier to leave it in) to enable switching

# RAF[, 'm', 'intervention', 1] = freqs # dim1 :n_insecticides # dim2 :sex: m,f # dim3 :site: intervention,refuge # dim4 : generations

see run_rot for initial test mortality calculation

Removed these comments from todo because think we have moved on : 1. check out these runs that have the mut-sel balance option enabled df_res_all_rot_cost10000_mut_sel_bal.rda

~ check on the runs with the freq limit removed : df_res_all_rot_cost10000_no_freq_limit.rda

~ think about whether I can implement allele frequencies just being constrained to their starting values the first time that each insecticide is used.

14/2/2019

starting to get mortality recording working currently looking at the df_mortali array from within a stop put inside run_rot()

with current starting resistance frequencies and effectiveness the mortalities start below the threshold where we would stop using the insecticide (90%)

we can set starting frequencies lower and effectiveness higher e.g. run_rot(eff=1, start_freqs=0.001)

moved rotations repo back from Ian

https://help.github.com/articles/about-repository-transfers/ Transfer is at the bottom of settings, options : from Ians account

20/2/2019

watched Ellies talk ~ she particularly talked about the benefit of combining complementary interventions in future (e.g. community larvel source management, ATSBs, etc.) ~ How these combine will depend on which part of the transmission cycle they target, e.g. outdoor biters ~ she showed quite nice cartoons of how untreated nets, treated nets and IRS can reduce tranmsission in different ways ~ made me think there may be potential to develop a clearer visualisation of the transmission cycle and how interventions can work, similar to stuff I did with Gerry ~ other interesting thing that Hilary asked about is how the different insecticides can modify both repellance and killing and how this may change due to resistance. Maybe similar to converstaions we had with Francis Rowland. ~ but I don't really have the inclination to move in the modelling direction but may be options to keep some interests going parallel with it.

~ more on how to record and display mortality

currently resistance frequencies : stored in an array (RAF), saved to a dataframe (df_res) converted to a diff format (df_res2) plotted with rot_plot_resistance()

for mortality : store in (a_mort) saved to df_mortali ...

Alternative Options A. convert df_mortali to df_mort2 so that can use rot_plot_resistance()

B. save resistance and mortality in same structure, maybe able to add mortality on as a column to df_res2

Tried the above, seems to work.

This shows there is a very consistent relationship between mortality & resistance. plot(df_res2$mortality ~ df_res2$resistance)

So maybe we don't even need to save mortality ???

this showed not quite linear plot(lm(df_res2$mortality ~ df_res2$resistance))

So can we just calculate a conversion from resistance to mortality that will allow us to estimate the resistance frequency that corresponds to 90% mortality ??

I can probably just calculate the conversion from the code I've written recently ?

Does this mean that the frequency differences of the different genotypes are not directly important ?

see genotype_freq() and mort_from_resist()

sent msg to Ian

21/2/2019

Cool getting there with mortality checking : Best to set logy to F to be able to see 90% mort

dfr <- run_rot(start_freqs = 0.0001, eff=1, logy=FALSE)

26/2/2019

starting to get mortality switch working.

seems OK from messages but plots don't look right

in the plots looks like still responding to frequency rather than mortality ??

Something a bit weird going on ... Currently the threshold mortality is not reached in the plots, but seems that frequency is.

Fixed!!! mostly a problem that eff,dom_sel,rr was not passed to new calc functions (and I didn't notice because they had defaults).

dfr <- run_rot(start_freqs = 0.000000001, eff=1, logy=FALSE, mort_or_freq='mort', rot_interval=0, threshold=0.9)

1/3/19

Had idea to try to plot regio of parameter space where mortality less than 90%. Not sure it is very useful, but did make a start.

see plot_mort_thresh.r

5/3/19

~ added mort_or_freq switch to rotresist1/ui.r ~ added threshold to UI

7/3/19 Using mortality switch for sequences, currently it switches back to other insecticides even when they are below threshold.

Aha I think that is due to the min_rwr_interval. When it is set to 0 then the sim doesn't proceed.

~ struggled to add legend for green mortality lines fixed now

8/3/19 eeek, found that there is a cran package called rotations, and this seems to be messing up deployment to shinyapps but was OK when I installed my rotations from github

~ published to shinyapps and sent link to Ian

12/03/19

~ change origin of git repo in RStudio, probably need to do from git remote: This repository moved. Please use the new location:
remote: git@github.com:AndySouth/rotations.git
from the commandline : git remote rm origin git remote add origin git@github.com:AndySouth/rotations.git

I think this is fixed now : Warning: Error in raf_get: object 'RAF' not found run_rot() can create it but when I add breakpoints it disappears

where is raf_get called ? insecticide_check() rot_migrate()

I suspect problem is here in insecticide_check()

end of rotation reached

if (gens_this_insecticide >= rot_interval) 
  {
  # 31/7/18 I could add a check in here for case
  # where only this insecticide remains below resistance threshold
  # in which case don't want to change
  # if I reset gens-this-insecticide then it will be used for another rot_interval
  #I can use either of these to get the freqs for all other insecticides
  # RAF[-current_insecticide, 'f','intervention',1]

e.g. raf_get(a) Error in raf_get(a) : object 'a' not found

It is only from scenarioA. Could it be something about sequence versus rotations ? Yes it doesn't happen when rot_interval=0 but did happen when rot_interval=10. But then when I start trying to look at it even from a new session it seems to dissappear.

14/06/2019

Testing before sending to Sam. Error in mortRR["f", "intervention"] : incorrect number of dimensions

This is when coverage set to 1, because there is no intervention dimension in mortRR[] in run_tot()

fixed

then got this from the script because it is trying to cast a tibble with an insecticide column and gens to a number Error: (list) object cannot be coerced to type 'double'

so I modified gens_under_thresh() to just return the numeric column

BUT this doesn't solve and I don't know why this problem appeared when it wasn't there before ???

gens_under_thresh() is giving one row per insecticide, is that what I'm expecting ?? No it is supposed to be : 'the sum across all insecticides for the number of generations below threshold while that insecticide is deployed'

so why has it changed ?? BUG introduced by the mortality addition

14/6/19 corrected bug here by adding back res <-

res <- summarise(res, tot_gens_dep_under50 = sum(.data$gens_dep_under50)) %>%

20/8/19

checking this message below from my notes, problem is that simulation continues even when mortality starts below thresholds. Inputs are set for sequence, what it does is switch between the insecticides for the min_rwr_interval. If the min_rwr_interval is set greater the length of the use periods increases too.

Probably I need to check that mortality not below threshold at start of the r_w_r interval ? In insecticide_check()

~~ why in this scenario does it rotate insecticides even though it starts below mortality threshold (to do with min_rwr threshold) run_rot( n_insecticides = 3 , max_gen = 150 , start_freqs = 1e-04 , rot_interval = 0 , eff = 0.9 , dom_sel = 0.5 , dom_cos = 0.5 , rr = 0.5 , expo_hi = 0.5 , coverage = 0.6 , migration = 0.01 , cost = 0.01 , logy = FALSE , no_r_below_start = FALSE , min_rwr_interval = 5 , threshold = 0.9 , mort_or_freq ='mort' )

this is the part in insecticide_check() that causes problem

if ( check_value > threshold &
     #20/8/19 also this allows an inseticide to start being used even when outside threshold
     gens_this_insecticide > min_rwr_interval )
{
  change_insecticide <- TRUE

OR does it ? surely the problem is that insecticide_switch() is choosing an insecticide that is already outside threshold ?

something weird going on in insecticide_switch() mortality threshold should be converted to survival for the comparison but doesn't seem to be :

switch from insecticide 1 to 3; new mortality vs threshold = 0.100044 vs 0.900000

AHA! NASTY BUG! there are multiple insecticides and it cycles through them because I do threshold <- 1-threshold that can give uncertain results.

4/9/19

However I change the run workflow I should try to keep it so that I can do the runs on multiple cores.

5/9/19

refactoring inputs Currently the list linputs which has one row per run is created in sensi_an_rotations1.Rmd by hardcoded min-max limits

new inputs plan Have a rotations function that accepts a csv file or defaults and uses that to create the linputs list, also it can save the csv. Closer to what resistance code did.

How to name the different files ? what are similar names in openMalaria ?

6/9/19

GOOD getting there, now a csv file saved in extdata that specifies the experiment

linputs_multi <- set_run_inputs()

to get all inputs for a single scenario

linputs <- purrr::map(linputs_multi, 2)

run one scenario

dfres <- do.call(run_rot, linputs)

10/9/19

got different insecticides working added get_one_in() to get inputs for one scenario from multi-senario file

11/09/19

This is the new workflow.

read an experiment file specifying ranges

inex <- read_in_expt() #specify filename

create an object containing inputs for all scenarios

linmulti <- set_run_inputs(inex=inex)

get inputs for a single scenario

linputs <- get_one_in(linmulti, scen_num=2)

run one scenario

dfres <- do.call(run_rot, linputs)

12/09/19

sorting mort_or_freq and other inputs in set_run_inputs() and get_one_in()

something I've done means that now the first insecticide is not changed in any scenarios...

linputs <- get_one_in(linmulti, scen_num=1) linputs <- get_one_in(linmulti, scen_num=2) linputs <- get_one_in(linmulti, scen_num=3)

linputs <- get_one_in(linmulti, scen_num=4) linputs <- get_one_in(linmulti, scen_num=5) linputs <- get_one_in(linmulti, scen_num=6)

run one scenario

dfres <- do.call(run_rot, linputs)

13/09/19

inex <- read_in_expt("_in_expt_rotations_default_one_scenario.csv", nscenarios=1) linmulti <- set_run_inputs(inex=inex) linputs <- get_one_in(linmulti, scen_num=1) dfres <- do.call(run_rot, linputs)

moving to have male_expo_prop default as 0

#now set to log-uniform 
minlog <- log(inex$start_freqs_min,10)
maxlog <- log(inex$start_freqs_max,10)

linmulti$start_freqs[i] <-  10^runif(1, min=minlog, max=maxlog)

Use hist to display the distribution of values.
hist(linmulti$start_freqs, breaks=100)

2/10/19

corrected bugs in rotation_runs2019.Rmd to do with saving rot_interval

3/10/19

rotation_runs2019.Rmd - working for 1000 runs but progress messages not appearing in terminal If I knit I do get progress messages and also get a document with all plots but not sure how that will work if running for 10k

~ improved plots so seq and rot on same page ~ changed gens_under_thresh to annotation to make text better

~ gens_under_thresh seemed to be counting 11 for a 10gen rotation ? Seems that a 10 gen rotation doesn't switch until generation 11 etc. ~ AND seemed that rot_interval not passing to run_rot ?

BOTH of these were because exit_rot default was true, now changed default to be false

TWO SEEMING PROBLEMS prob1) sequence insecticide continues to be used after threshold I expect may be a problem with exit_rot=FALSE prob2) migration in the new runs seems to make active & refugia the same, maybe I've set limits too high ? I have default migration set to 0.01 and limits currently set to 0.1 to 0.9 df_res2 <- run_rot(rot_interval = 0, mort_or_freq = 'freq', migration=0.5)

Having difficulty replicating prob1

df_res2 <- run_rot(rot_interval = 0, mort_or_freq = 'freq', migration=0.5)

I set min_rwr_interval to 10 to stop rapid switches for sequence.

run_rot(n_insecticides=2, mort_or_freq='freq', cost=0.013, start_freqs=0.077, rot_interval=10, eff=0.98, dom_sel=0.55, dom_cos=0.55, rr=0.47, expo_hi=0.69, coverage=0.23, migration=0.86, max_gen=500, min_rwr_interval=1)

note there is a problem with the quotes when copying from the pdf

diagnostics not helpful

Aha, think its just to do with exit_rot = FALSE run_rot(n_insecticides=2, mort_or_freq='freq', rot_interval=10)

4/10/19

Fixed rotation insecticide continues to be used after threshold which was a bug in insecticide_check() introduced when mort thresholds added

NEW SEEMING PROBLEM run_rot(n_insecticides=2, mort_or_freq='freq', rot_interval=10, exit_rot=FALSE, diagnostics = TRUE, min_gens_switch_back = 0) Aha I think problem that insecticide itself is not assessed so keeps being used even when abovbe thresh FIXED

fixed gens_under_thresh count which had changed to be for all insecticides, not just one.

successfully GOT runs2019 to access the xls file for experiments directly

7/10/19

got diff insecticide runs working Got mortality to appear properly in plots & key in non dispersal runs.

8/10/19

added to readme

9/10/19

spent the day sweating over getting migration working.

Now I'm pleased I think it is simpler and I know what its doing. Migration 1 is random mixing.

see the excel file for testing : rotations-dispersal-demo-andy201910.xlsx

10/10/9

migration effects on gens under thresh run_rot(coverage=0.5, cost=0, migration=0) #64 run_rot(coverage=0.5, cost=0, migration=0.01) #72 run_rot(coverage=0.5, cost=0, migration=0.1) #92 run_rot(coverage=0.5, cost=0, migration=0.3) #104 run_rot(coverage=0.5, cost=0, migration=0.6) #116 run_rot(coverage=0.5, cost=0, migration=0.8) #120 run_rot(coverage=0.5, cost=0, migration=1) #124

migration limits ? 0.01 - 0.6

But other A3 scenarios showe that nearly all runs with migration > 0.2 or 0.3 have treated area & refugia with very similar frequency curves.

I noticed that male_expo_prop min was 1 in earlier runs, changed to 0

14/10/19

Quotes from IRM in Aedes paper : Dusfour2019 Management of insecticide resistance in the major Aedes vectors of arboviruses: Advances and challenges.

in blue could be disagreed with :

Box 1. IRM strategies Rotation

This strategy employs temporal alternation of insecticides with different modes of action and is based on the assumption that resistance diminishes insect fitness in the absence of insecticide.

[in our rotation paper we could address this misconception that fitness costs are required for rotations to be beneficial].

All comparisons with mosaic, rotation, or mixture based on theoretical models or empirical studies demonstrate that “reactive alternation” underperforms the other IRM approaches in slowing down resistance evolution because alternation is implemented too late, when resistance alleles are already at high frequency [60].

In Singapore, permethrin was replaced by pirimiphos-methyl for thermal fogging after reduced efficacy was observed [61], yet after 9 years of interruption, high permethrin resistance persisted in field populations.

Persistence of high levels of pyrethroid resistance was also observed 7 to 10 years after they were replaced by organophosphates in the state of Sao Paulo (Brazil) [62]. As mentioned before, however, maintenance of resistance can also be due to the exposure of vector populations to household and agricultural insecticides or other xenobiotics (pollutants, other pesticides, etc.) [63–65].

On the bright side, the Onchocerciasis Control Programme in West Africa is probably the most emblematic case of successful IRM in public health. After the first detections of temephos resistance in Simulium damnosum, a rotation scheme of insecticides with alternative modes of action was implemented for larval treatments. The rotation strategy drove the reversion of organophosphate resistance and preserved S. damnosum insecticide susceptibility for more than 20 years until the program was stopped [66].

For malaria, a large-scale trial was conducted in Mexico to evaluate the impact of different IRM strategies on the resistance of Anopheles albimamus [67]. In this study, a 3-compound annual rotation, a 2-compound mosaic, and a long-term treatment with a single insecticide were compared. An increase of resistance was observed in all villages, but it demonstrated that mosaic or rotation selected for low resistance levels and remained stable compared with villages with single treatment where resistance increased rapidly and to levels significantly greater [68]. Control failures associated with pyrethroid resistance in malaria vectors have been observed for IRSs using pyrethroids. Switching to DDT against An. funestus populations without kdr mutation but specific metabolic mechanisms of pyrethroid resistance in South Africa [69] or to bendiocarb against An. gambiae in Equatorial Guinea restored the control of resistant populations [70]. ... The limited evidence-based information to support any IRM strategy in public health (e.g., mixture, rotation, noninsecticidal tools, etc.) complicates the prioritization of interventions.

Despite that, we have to advocate for the introduction of IRM policies into vector control programs by favoring a proactive approach over a response mode.

For example, it has been demonstrated that the preset rotation of unrelated insecticides performed better than changing insecticides over time in response to resistance [60, 66].

[next sentence is good for justifying future work]

Until now, it has been near impossible to anticipate how resistance will impact the efficacy of vector control tools in the field and how a given IRM strategy will slow or reverse the evolution of resistance in targeted populations. A continuous effort connecting research institutes, control programs, industry, and WHO is required to conduct operational research that evaluates vector control strategies under field conditions and the genetics of resistance in broad terms (genotype–phenotype relationships, mosquito behavior, modeling). The outcomes of these studies will provide crucial information to optimize decision-making systems and the procedures to adjust IRM strategy according to field situations.

15/10/2019

seeming bug with mortality switch that for rotations it stops the run when one insecticide is still below threshold.

e.g D4 scenario1 run_rot(rot_interval=10, n_insecticides=2, start_freqs=0.0018, cost=c(0.083,0.047), eff=c(0.95,0.96), dom_sel=c(0.18,0.41), dom_cos=c(0.85,0.98), rr=c(0.29,0.71), expo_hi=0.69, male_expo_prop=0.17, coverage=0.86, migration=0.28, min_rwr_interval=10, mort_or_freq = 'mort')

but compare with freq that does work :

run_rot(rot_interval=10, n_insecticides=2, start_freqs=0.0018, cost=c(0.083,0.047), eff=c(0.95,0.96), dom_sel=c(0.18,0.41), dom_cos=c(0.85,0.98), rr=c(0.29,0.71), expo_hi=0.69, male_expo_prop=0.17, coverage=0.86, migration=0.28, min_rwr_interval=10)

in insecticide_check()

I think I fixed this issue ...?

6/11/2019

starting to get some of the analysis code ready for Sam.

taking code from sensi_an_rotations1.Rmd

13/11/2019

Why does this sequence switch back when there is no cost or dispersal ?

run_rot(rot_interval=0, n_insecticides=3, start_freqs=0.0021, cost=c(0), eff=c(0.98), dom_sel=c(0.51),dom_cos=c(0.63), rr=c(0.8), expo_hi=0.47, male_expo_prop=0.01, coverage=1, migration=0,min_rwr_interval=10, diagnostics = FALSE)

Using a mortality switch it now doesn't switch back : run_rot(rot_interval=0, n_insecticides=3, start_freqs=0.0021, cost=c(0), eff=c(0.98), dom_sel=c(0.51),dom_cos=c(0.63), rr=c(0.8), expo_hi=0.47, male_expo_prop=0.01, coverage=1, migration=0, min_rwr_interval=10, mort_or_freq = 'mort')

Not due to dom_cos >0 even though cost=0.

insecticide1 check=0.488229891221731 freq=0.488229891221731 surv=0.406690062481814 gens_this_insecticide=39 generation 41

insecticide1 check=0.520515732451062 freq=0.520515732451062 surv=0.431997734603674 gens_this_insecticide=40 change from insecticide1 check=0.520515732451062

generation 121, switch from insecticide 3 to 1; new frequency vs threshold = 0.488603 and 0.500000

Seems like freq has declined from 0.5205 to 0.4886 while insecticide not in use despite no costs or dispersal ???

ooo weird seems to go down the generation after use stopped and then stays at that constant level:

gen39 i1 freq=0.455590826478509 gen40 i1 freq=0.488229891221731 gen41 i1 freq=0.520515732451062 gen42 i1 freq=0.488603067950344 gen43 i1 freq=0.488603067950344 gen44 i1 freq=0.488603067950344

Ian suggested this could be because we assess frequency after selection and before random mating. So what could be happening is the sequence is stopped because f is above threshold, then random mating happens to take it below threshold.

In which case we could just make sure that insecticide_check is called after mating. But was tricky to do, partly a result of us assessing just Females for insecticide switch. Now the switch is checked based on mean of M&F.

New setting to make it easier to bug-check

run_rot(rot_interval=0, n_insecticides=3, start_freqs=0.1, max_gen = 50, cost=c(0), eff=c(0.98), dom_sel=c(0.51),dom_cos=c(0.63), rr=c(0.8), expo_hi=0.47, male_expo_prop=0.01, coverage=1, migration=0,min_rwr_interval=10, diagnostics = FALSE)

*storedgen16 freqs i1:0.5172 i2:0.1 gen17 freqs i1:0.4852 i2:0.1292 gen17 i1 freq=0.485228904909616

19/11/2019

confirm that insecticide check based on mean of M&F does what I'd expect

trying 100 runs of A1

NO sequence still not doing what I' expect, i.e. for A1 it's going back to insecticide 1 : Ah probably just because I haven't built and saved ...

run_rot(rot_interval=0, n_insecticides=2, start_freqs=0.077, cost=c(0), eff=c(0.98), dom_sel=c(0.83), dom_cos=c(0.47), rr=c(0.2), expo_hi=0.69, male_expo_prop=0.17, coverage=1, migration=0, min_rwr_interval=10)

Currently I have min_gens_switch_back and min_rwr_interval set to 10. Probably I should set them to 0.

Instead have an input specifiying interval on which insecticides can be changed - applying both to sequence and rotation.

change_interval = 10

if (gens_this_insecticide %% change_interval == 0)

FIXED and new 16k simulations run.

25/11/19

paper from Justin :

Busi 2019 Rotations and mixtures of soil-applied herbicides delay resistance

From the abstract : Simulations indicate that without rotation it takes twice as long to select for resistance to a particular soil-applied herbicide – trifluralin – than to any other herbicide option considered. Relative to trifluralin-only use, simple herbicide rotation patterns have no effect in delaying resistance, whereas more complex rotation patterns can delay resistance two- or three-fold. Herbicide mixtures further delay resistance up to six-fold in comparison to single use or simple herbicide rotations.

The model is able to represent the evolution of resistance conferred by either single or multiple genes (monogenic or polygenic resistance) using explicit individual-based genetics

The PERTH model, parametrized with data obtained from a decade of experimental evolution studies,14,30–32,35,42 predicts that herbicide mixtures are the most robust strategy to minimize the risk of resistance and cross-resistance evolution to soil-applied herbicides. The relative benefits of rotations and particularly mixtures in delaying resistance were found to be consistent across a number of genetic and fitness-penalty scenarios and 100-fold difference in initial allele frequencies.

Unsurprisingly, the worst ‘strategies’ were those with no herbicide rotation (1111 best, 2222, 3333).

Table 2. The rank column indicates the relative efficacy of the strategy in delaying resistance under GS1, which is not affected by allele frequency GS3 versus GS1.

Pesticide mixtures at high dosages of each pesticide and/or combinations appears to be the most effective measure at delaying resistance evolution because they help ensure that pests are treated with two different pesticides. Similarly, weeds resistant to one herbicide mode of action are likely to be killed by a another herbicide acting at a different site of action.

Me : (not sent)

Hi Justin,

Thanks for sending this on, apologies it has taken me a whole month to get around to reading it in detail. I've just submitted th final proofs for the windows of selection paper and have finalised the simulations for our rotations paper.

It was encouraging to see some plots that looked similar to ours.

However I do find this paper somewhat confusing. The title 'Rotations and mixtures of soil-applied herbicides delay resistance' suggests a benefit to rotations.

Then in the abstract it says : "Relative to trifluralin-only use, simple herbicide rotation patterns have no effect in delaying resistance, whereas more complex rotation patterns can delay resistance two- or three-fold."

I wasn't able to establish what the mechanistic difference was between the simple and complex rotations.

They do seem to show results and justification for mixtures being beneficial.

I think, looking at Table 2, that they are comparing rotations to constant use of a single herbicide. I don't think that they consider sequential stratgies (i.e. switching insecticides once resistance thresholds are reached).

In our work, that we are close to submitting. We find very little difference between sequential and rotation strategies.

Andy

So for one 'good' herbicide, simple rotations provide no advantage. Complex rotations apparently can provide an advantage, but I wasn't able to establish why from their results. They results do seem to provide clear support and justification for the benefit of mixtures.

27/11/2019

paper_figs_rotations201911.Rmd

choose example runs for early figs

ARg, this may not matter, but I had wondered if I could directly compare runs with and without costs. BUT it seems the way that random samples are drawn that th values for other input change.

e.g. on p6, A1 no cost eff=0.79 A2 with cost cost=0.79, eff=0.99

run_rot(rot_interval=0, n_insecticides=5, start_freqs=0.092, cost=c(0), eff=c(0.79), dom_sel=c(0.16), dom_cos=c(0.94), rr=c(0.19), expo_hi=0.73, male_expo_prop=0.15, coverage=1, migration=0, min_rwr_interval=10)

run_rot(rot_interval=0, n_insecticides=5, start_freqs=0.077, cost=c(0.079), eff=c(0.99), dom_sel=c(0.5), dom_cos=c(0.81), rr=c(0.38), expo_hi=0.46, male_expo_prop=0.16, coverage=1, migration=0, change_interval=10)

This may not matter for the main runs. but may mean that I want to choose example runs differently.

Happens in set_run_inputs()

# select random numbers for each insecticide
cost <-     runif(n_rands, min=inex$cost_min,        max=inex$cost_max) 
eff  <-     runif(n_rands, min=inex$eff_min,         max=inex$eff_max)

Interestingly, (and possibly worryingly) seems that the random deviate isn't used if min and max are the same.

9/12/19

Thinking that I can output mean % mortality over all generations and that will give us a better comparison between sequences and rotations (alongside going back to recording total generations untill all insecticides exhausted).

10/12/19

corrected recordning of gens_under_thresh for mortality even though probably not going to use anymore

~ set recording of i total generations (maybe from run_rot) ii cumulative or mean mortality (maybe post-processed)

I think better to use mean mortality because it is more obviously meaningfull, and a potential way of comparing runs even if they have different lengths or number of insecticides.

2020-01-06 NEW YEAR

sorted 3 outputs : gens_rot_minus_seq #total generations gens_u_rot_minus_seq #gens under threshold mort_rot_minus_seq #average mortality

problem that gens_u_rot and gens_u_seq both currently 0

seems to be because linputs$freq_thresh & mort thresh are both NULL

This B2 scenario 33 p34 is interesting, big benefit of rotations. One insecticide stays low in rotation seemingly due to use of others.

run_rot(rot_interval=0, n_insecticides=3, start_freqs=0.013, cost=c(0.061,0.024,0.047), eff=c(0.97,0.88,0.72), dom_sel=c(0.49,0.15,0.92), dom_cos=c(0.72,0.8,0.37), rr=c(0.25,0.51,0.42), expo_hi=0.52, male_expo_prop=0.92, coverage=1, migration=0, max_gen=500, change_interval=10)

saved run file for sam to rotations_runs2020_sam.rmd

2020-01-10

sam did new runs

in paper_figs_rotations201911.rmd ~ replace rot_minus_seq with gens_rot_minus_seq

2020-01-13

looking at mean mortality from new runs

2020-01-17

A2 runs run 94 (p95) rotation better p109 sequence better

looking for the B1 runs where rotations are a bit better than sequences

2020-02-25

checking on B1 runs following feedback from Sam.

B1 scenarios, (no costs, no dispersal, insecticides different).

In earlier runs there was no advantage to rotation when there is no cost and dispersal. In the latest runs there is.

I think this is something to do with time spent above the resistance threshold. You may remember that insecticides can only be changed on a yearly timestep. This can result in resistance going above threshold, when that happens it can result in simulations lasting longer.

It seems that something is happening with this that gives the advantage to rotations. This creates a slightly strange 'dog-tooth' pattern seen in fig x3 and pasted below.

But sam checked that in the example 200 runs there were just a few that had minor benefits for sequences (not rotations)

I could try doing 1000 runs to check ? or do an equiv of fig x3 for the 200 runs.

its x2 which has the dogtooth OR first check how the fig x2 is produced to check for possible glitches

fig x1 shows that of the runs that are different they mostly show benefit of rotation.

fig-si-1-abs-diff-rot-seq shows NO B1 runs with benefit of sequence.

2020-02-25 checking B1 scenarios that are giving weird results

df_res_B1 <- dplyr::filter(df_res_16, idAD=='B' & id14 == 1)

df_res_B1 %>% group_by(gens_rot_minus_seq) %>% tally()

A tibble: 7 x 2

gens_rot_minus_seq n

1 -10 1

2 0 6454

3 10 1653

4 20 1207

5 30 514

6 40 150

7 50 21

so this seems to suggest that rotations are indeed better for B1

possibilities :

A something different in the runs that Sam did and the 200 example ones done for the document

B a problem with the document

but the B1 document looks as expected ???

should run 1,2,3 etc. from the doc be expected to be the same as the ones in df_res_B1 ?

remember pdfs with figs are in C:\Dropbox\resistanceResults\rotations

YES for run 1 the inputs look about the same

EXCEPT that e.g. eff is recorded in df_res_B1 as just a single value (0.98) whereas there are 2 insecticides should be c(0.98,0.74)

also starngely the pdf seems still to be showing gens in-use under thresh (maybe I didn't get around to changing that)

a) install-restart and try redoing B1 pdf for first 10 runs 1) check on changing gens-under-thresh output for pdf 2) check on saving multiple values when insecticides are diff

run1 in pdf rot:40, seq:40 in runs rot:50, seq:40

I wonder if something went wrong with Sams runs that just running 1 insecticide (but then would expect them to be just like A1 runs ??)

& maybe Sams reporting on the pdf was based on the gens-under-thresh numbers which are now superceded (my mistake for not updating package before re-running pdf)

Can I run the rot for run1 from the pdf, does it generate 40 as in the pdf or 50 as in the runs ? eek maybe something hadn't updated before Sam did the runs ?

run_rot(rot_interval=10, n_insecticides=2, start_freqs=0.077, cost=c(0,0), eff=c(0.98,0.74), dom_sel=c(0.55,0.55), dom_cos=c(0.24,0.76), rr=c(0.77,0.47), expo_hi=0.69, male_expo_prop=0.17, coverage=1, migration=0, max_gen=500, change_interval=10)

this gives 40 gens, same as the pdf.

Now can I rerun the first 10 runs in the same way that Sam did ? To do that can just click on the saved rda file. It does give the same results as the pdf.

I think I may need to ask Sam to re-run. Making sure that we've updated everything.

One remaining query, why are the inputs for multiple insecticides not saved into the df_res_all object (by both mine & Sams runs) ?

it must happen here in rotation_runs2019 : df_in_out <- as.data.frame(c(linputs, gens_rot=gens_rot, gens_seq=gens_seq, gens_u_rot=gens_u_rot, gens_u_seq=gens_u_seq, mort_rot=mort_rot, mort_seq=mort_seq, id=i))

linputs has the multiple entries for eff, df_in_out has lost them i thought I had done something to allow for that in past ?

(to be clear I don't think the above is source of an error, it just might influence ability to analyse runs with diff insecticides later)

BUT e.g. when there are 4 insecticides df_in_out ends up with 4 rows ... & the outputs repeated THEN this bit (i think) just saves one row so the outputs will be OK

df_res_all[i,] <- df_in_out

It's a bit of a mess though.

#2020-02-25 BEWARE (see above) when insecticides are different this saves just first row of df_in_out #result is that correct results are saved but with just the inputs for the first insecticide #this shouldn't effect the output results, but I should try to correct #putting results into one row of overall results dataframe df_res_all[i,] <- df_in_out

this demonstrates problem

tst <- as.data.frame(c(linputs))

How did I solve it somewhere else ?

solution may be to use linmulti which contains the strings with multiple insecticide values in. e.g. these are the effectivenesses for the first run linmulti$eff[1] [1] "0.98304248759523, 0.738747693016194"

get_one_in() converts linmulti to linputs which can be used to run a simulation, but maybe I just need to take a single row from linmulti to save into the output file.

Maybe add an argument to get_one_in to be able to keep the strings ?

tst <- as.data.frame(c(linputs_tosave))

but still remaining problem that df_res_all can't store a string so it seems to store 1 instead, e.g. for eff

2020-02-25 : get Sam to rerun the 16k runs

2020-10-05 reviewing runs for rotations paper & Qs from Ian

Andy….I assume the following description is how you transformed the data to plot Figure 7 ? We quantify the relative performance of sequential and rotations policy (see later discussion of Figures 7 and 8) using the metric, z. z=(L(r)/L(s) -1)*100

this is my fig x1, the output is rms_percent

df_res_16$rms_percent <- 100 * df_res_16$gens_rot_minus_seq / df_res_16$gens_seq

Ian asks to change the mortality fig (8) to %difference.

I think it is : figxmort-which-best-meanmort

The rational for keeping it to absolute mortality differences is that the values are meaningful to those interested in insecticide resistance, and that mortality is already a %, so talking about percent of a percent gets really confusing. Currently the figure shows that most runs result in just a 1% difference in mortality - that readers may equate e.g. to the difference between 89 & 90% in a bioassay.

TODO 2020-20

TODO (from before)

Try to fix so that this saves multiple insecticide inputs, e.g. as c(x,y,z) df_in_out <- as.data.frame(c(linputs, gens_rot=gens_rot, gens_seq=gens_seq, gens_u_rot=gens_u_rot, gens_u_seq=gens_u_seq, mort_rot=mort_rot, mort_seq=mort_seq, id=i))

paper_figs_rotations201911.Rmd see rotation_runs2019.Rmd & rotations-experiments-creating-201909.xlsx

convert single run plots to output total generations rather than gens-under-thresh

check on dispersal limits using new method. Seems that anything above 0.4 leads to near complete mixing run_rot(migration=0.4, coverage=0.5) interesting how small levels of migration cause resistance freqs to decline even with no costs write word description of what dispersal does rename migration to dispersal ? migration limits ??? I need to check on them more ! to make sure I know what they are doing.

~ consider renaming package to avoid conflict with CRAN rotations package insectrot rotinsect rotseq rotateinsect

CHECK ALL INPUT PARAMS

RUNS FOR PAPER September 2019

A identical insecticides, seq vs rot, resistance frequency switch, A1 no costs no dispersal, A2 costs no dispersal, A3 dispersal no costs, A4 costs & dispersal. Sensitivity analysis of inputs.

B different insecticides, as above, no sensitivity analysis.

C identical insecticides, mortality switch, reduced input range so mortality starts above 90%

D different insecticides, mortality switch

Runs are combinations of : Strategy (2) : sequence, rotation Insecticides (2) : same, different Switch (2) : resistance frequency, mortality Scenarios (4) : no costs no dispersal, costs no dispersal, dispersal no costs, costs & dispersal Replicates (10000) :

So number of runs 222410000 = 320000

STill if have min_rwr_interval > 0 an insecticide at very start of simulation outside the threshold will continue to be used. Mostly I think min_rwr_interval not a good idea, but does kind of represent what you might expect to happen in field when it is set to a year i.e. that once sprayed can't be unsprayed.

dfr <- run_rot(start_freqs = 0.000000001, eff=1, logy=FALSE, mort_or_freq='mort', rot_interval=0, threshold=0.9)

~ fix that mortality calc in plots and that for switching is currently repeated ~ check it works with diff insecticides i.e. eff=c(0.7,0.8,0.9) I suspect not because mort_from_resist() or something else not properly vectorised.

protocol for sensitivity analysis ~ store input ranges ~ generate inputs for each run ~ run ~ store outputs ~ calculate output summary stats per run ~ concatenate inputs & outputs ~ relate outputs to inputs

~ check in resistance::fitnessSingleLocus() on separation of dominance of cost & selection. Was probably OK for Levick et al because we didn't have cost in. But I think that dominance of cost may be set to 0 which means that the RS experiences no effect of cost.

a_dom       <- array_named(locusNum=c(1,2), exposure=c('no','lo','hi')) 
a_dom[1, 'hi'] <- dom1
a_dom[2, 'hi'] <- dom2

#exposure 0 'no'
a_fitloc[ paste0('RS',locusNum), 'no'] <- 1 - (a_dom[locusNum, 'no'] * a_cost[locusNum])
a_fitloc[ paste0('RR',locusNum), 'no'] <- 1 - a_cost[locusNum]

state of play

~ sensi_an_rotations1.Rmd ~~ produces an output dataframe containing both inputs and outputs ~~ saves outputs to df_res_all_rot.rda

~ rotations1.Rmd : example plots ~ UI : rotations/shiny/rotresist1/server.r https://andysouth.shinyapps.io/rotresist1/ ~ RAF[,,,2] to get allele frequencies for generation 2 (by insecticide,sex and active|refuge)

useful

~ one generation per month, 12 (or 10) generations per year, 500 gens is ~40-50 years

ian-hastings/rotations documentation built on Dec. 14, 2020, 11:42 p.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.

ian-hastings/rotations Simulating the Evolution of Insecticide Resistance under Rotations

In ian-hastings/rotations: Simulating the Evolution of Insecticide Resistance under Rotations

rotations_notes.Rmd

' n_insecticides <- 3

' # list(c()) is important below to create one list for each dim

' l_dimnames <- rep(list(c('no','lo','hi')), n_insecticides)

' names(l_dimnames) <- paste0("niche", 1:n_insecticides)

' do.call(array_named,l_dimnames)

this feels hacky, but it does work & I can put in a function

create the flexible r* columns

create inflexible columns ready to add r* columns on

diff frequency curves for r1 & r2

high cost, no migration

no cost, no migration

no cost and rotate when resistant, no migration

no cost, some migration

testing multiplying RAF arrays

to allow making code more concise

programing with dplyr

resistance files to do runs

Q for Ian :

discussion

rotation :

sequence :

scenario 79 expo_hi 0.79 eff 0.72 rot_interval 0

tot gens deployed under freq 0.5 = 199

scenario 79 expo_hi 0.79 eff 0.72 rot_interval 6

tot gens deployed under freq 0.5 = 251

weirdly when I went to 3 insecticides, & cost 0.1

gen calc gives same result for each at 69 and plots finish at same point ?

aha! that was because max_gen=70 so max result could be 69 i.e. no insectidides in either scenario reached threshold

add in_use columns

remove NA refuge rows

linputs$migration <- 0.1

scenario 137 sum generations deployed insecticide under thresh rotation: 500 sequence: 258

scenario 143 sum generations deployed insecticide under thresh rotation: 500 sequence: 126

advantage to rotation despite dom_cos=0.23, cost=0.02

scenario 153 sum generations deployed insecticide under thresh rotation: 406 sequence: 289

30/4/18

setting params for 2167

this shows a small tweek of the params leads to the sequence lasting forever too

and I don't know which input change causes the seq to last forever ?

Aha seems that I hadn't included rr=0.78, rr default is 0.5,

with lower rr less selection for r and sequence does better

good refs for rotations paper

13/11/2018

23/11/18

5/12/18

6/12/18

6/2/19

13/2/2019

14/2/2019

moved rotations repo back from Ian

20/2/2019

21/2/2019

26/2/2019

1/3/19

5/3/19

12/03/19

end of rotation reached

14/06/2019

14/6/19 corrected bug here by adding back res <-

20/8/19

this is the part in insecticide_check() that causes problem

4/9/19

5/9/19

6/9/19

to get all inputs for a single scenario

run one scenario

10/9/19

11/09/19

read an experiment file specifying ranges

create an object containing inputs for all scenarios

get inputs for a single scenario

run one scenario

12/09/19

run one scenario

13/09/19

2/10/19

3/10/19

ian-hastings/rotations
Simulating the Evolution of Insecticide Resistance under Rotations