Introduction

Insect resistance to pesticides is a wicked and widespread problem (Gould, Brown, and Kuzma 2018). Predicting, identifying, and ultimately delaying the evolution of insecticide resistance is therefore critical for ensuring sustainable global food security (Hawkins et al. 2019; Lykogianni et al. 2021). To manage pesticide resistance, various strategies have been proposed to disrupt or weaken selection. These strategies often focus on varying the location and timing of pesticide application, and therefore varying the strength of selection for pesticide resistance, so that pest susceptibility to one or more pesticides is maintained (Ives and Andow 2002; Liang et al. 2013; Carrière, Fabrick, and Tabashnik 2016; Liang, Tang, and Cheke 2018; South and Hastings 2018; Sudo et al. 2018; Saikai, Hurley, and Mitchell 2021). Pesticide resistance management has mostly been associated with the effect of single genes. In such cases, resistance alleles have binary effects on phenotype, enabling resistant phenotypes to arise from genetic changes at single or a small number of loci (Raymond et al. 2001; Richard H. Ffrench-Constant, Daborn, and Le Goff 2004; Richard H. Ffrench-Constant 2013; Mangan et al. 2023). Nevertheless, resistance attributable to polygenic effects is also well-established (Via 1986; Firko and Hayes 1990; Busi, Neve, and Powles 2013; Haridas and Tenhumberg 2018; Green, Stenberg, and Lankinen 2020). And the relevance of polygenic resistance is likely to further increase given a rising interest in biological controls for sustainable crop protection (e.g., Green, Stenberg, and Lankinen 2020; Elnahal et al. 2022; Mangan et al. 2023). Developing strategies to maximise the efficacy of these tools is critical, and such strategies should be well-informed by predictions made from detailed, quantitative genetic models. Here we introduce the resevol R package as a tool for building individual-based models and simulating pest management (Daly et al. 2022).

The resevol package applies individual-based modelling and a quantitative genetics approach to simulate the evolution of a pest population on a changing landscape. Multiple traits determine the overall fitness of any pest genotype. For example, while alleles conferring resistance to a particular pesticide enhance fitness in the presence of that pesticide, they could also impose reproductive fitness costs. Such trade-offs can be quantified by the genetic covariance between traits. A focal goal of the resevol package is to model traits with a pre-specified, but potentially evolving, genetic covariance structure. To achieve this goal, each individual has a genome with L loci that underlie a set of T potentially evolving traits. Pleiotropic loci can vary in their effects on the direction and magnitude of polygenic traits, causing population-wide trait covariance to arise mechanistically from the underlying genetic architecture of individuals. To achieve this, two separate steps are necessary. First, an evolutionary algorithm is used to find a network of internal nodes that map standard random normal loci values to covarying traits. Values used to map loci to traits are incorporated into individual genomes. Second, a population of asexual or sexual individuals is initialised and simulated on a spatially explicit landscape separated into distinct units (e.g., farms). Land units can apply one of up to 10 pesticides, and one of up to 10 landscape types at a given time (for simplicity, here we interpret landscape types to be crop types). Pesticides and crops rotate independently within farms over time in a pre-specified way. The resevol package can thereby model complex and evolving agricultural pest traits over realistic landscapes that undergo different pesticide use and crop regimes.

Covarying pest quantitative traits

The first step of simulation is building individual genomes. This step is separate because it is computationally intensive, and genomes that are built might need to be inspected and stored. High computation time is due to the mechanistic nature of how genomes and covarying traits are modelled. Instead of imposing a trait covariance structure directly, an evolutionary algorithm is used to find a network that maps standard random normal values (loci) to covarying values (traits). This is useful because it allows genomes to model potentially evolving physiological constraints, and trade-offs among traits, from the bottom up. Since multiple networks can potentially map loci to the same trait covariance structure, it is possible to replicate evolution with different randomly generated genetic architectures. This approach to modelling individual genomes and traits thereby increases the complexity of questions that can be addressed for simulating evolution in agricultural pests.

Figure 1. Example network mapping loci (green circles) to traits (red diamonds) through an intermediate set of hidden layers (blue squares) in the mine_gmatrix function. Individual genomes in the resevol R package consist of standard random normal values for loci, real values for black arrows linking nodes, and real values for traits. Values shown for loci and arrows are an example for illustration.

Conceptually, the relationship between individual genotypes and traits is defined by a network connecting loci to traits through a set of hidden internal nodes (Figure 1). Values for loci are randomly drawn from a standard normal distribution, 𝒩(0, 1). Links between loci, internal nodes, and traits, can take any real value and are represented by black arrows in Figure 1. Traits can take any real value, and are calculated as the summed effects of all preceding nodes (i.e., the blue squares immediately to the left of the traits in Figure 1). Mathematically, loci are represented by a row vector of length L. Effects of loci on the first layer of internal nodes (black arrows emanating from loci in Figure 1) are represented by an L × T matrix, and transitions between internal nodes, and between the last set of internal nodes and the final traits, are represented by T × T matrices. For the example in Figure 1,

$$ \begin{pmatrix} 0.1, & -0.1, & 0.2 \\ \end{pmatrix} \begin{pmatrix} 0.4, & 0.1 \\ -1.2, & -0.6 \\ 0.2, & 2.2 \\ \end{pmatrix} \begin{pmatrix} 1.3, & -0.5 \\ -0.1 & 0.0 \\ \end{pmatrix} \begin{pmatrix} 2.2, & 0.9 \\ 1.6 & 0.8 \\ \end{pmatrix} = \begin{pmatrix} 0.2998, & 0.1081 \\ \end{pmatrix}. $$

Values mapping loci to traits become part of an individual’s genome, so the genome for the individual represented by Figure 1 is stored in the model as shown below:

0.1, -0.1, -0.2, 0.4, 0.1, -1.2, -0.6, 0.2, 2.2, 1.3, -0.5, -0.1, 0.0, 2.2, 0.9, 1.6, 0.8

Individuals with different loci can therefore have different covarying traits that are constrained by the network structure encoded in each genome. Individuals can be haploid (as in Figure 1) or diploid (in which case, allele values are summed at homologous loci).

An evolutionary algorithm is used to find appropriate values that produce covarying traits from loci (see Supporting Information 1 for details). Evolutionary algorithms are heuristic tools that can simulate adaptive evolution to find solutions for a broad range of problems (Hamblin 2013; Duthie et al. 2018). In the resevol package, the mine_gmatrix function runs an evolutionary algorithm and requires the argument gmatrix, which specifies the desired trait covariance matrix. The function initialises a population of npsize separate, and potentially unique, networks (i.e., npsize copies of a network like the one shown in Figure 1), and this population evolves until some maximum iteration (max_gen) or minimum expected network stress (term_cri) is met. In a single iteration of the algorithm, values mutate and crossover occurs between networks (i.e., some of the npsize networks swap values with some probability). Next, trait covariances produced for each network are estimated by initalising indivs individuals with loci sampled from a standard normal distribution. Network stress is calculated as the logged mean squared deviation between estimated covariances and those in gmatrix. Tournament selection (Hamblin 2013) is then used to determine the networks for the next iteration of the algorithm. Networks with the lowest estimated stress have the highest fitness, so these networks are disproportionately represented in the next iteration. Throughout the evolutionary algorithm, the lowest stress network is saved and returned upon network termination. The robustness of this network’s stress to sets of individuals with different loci values can be tested using the stress_test function.

An example run of mine_gmatrix with the same number of loci and internal nodes as in Figure 1 is shown below for traits that have an intended covariance of -0.4:

trait_covs  <- matrix(data = c(1, -0.4, -0.4, 1), nrow = 2, ncol = 2);
new_network <- mine_gmatrix(loci = 3, layers = 2, gmatrix = trait_covs, 
                            max_gen = 1000, term_cri = -6.0, prnt_out = FALSE);

The code above found a genome that produced the following expected trait covariances:

##            [,1]       [,2]
## [1,]  1.0008181 -0.3925799
## [2,] -0.3925799  1.0271340

The mean deviation between elements of the above matrix and the identity matrix provided by trait_covs is 2.1176024^{-4}. Lower values of term_cri and higher values of max_gen will result in a lower stress, but this will require more computation time, especially if the number of traits is high. Similarly, higher values of indivs will result in more accurate estimations of true stress, but this also requires more computation time. Additional arguments to mine_gmatrix can also be used to improve the performance of the evolutionary algorithm (see Supporting Information 1).

Simulating landscape-level pesticide resistance

The full output of mine_gmatrix is passed to the run_farm_sim function, which initialises and simulates an evolving population of pests on a changing landscape for any natural number of time steps (time_steps). In this section, we explain the landscape, pest ecology, and evolving pest traits.

Pest ecology

Individual pests can be modelled to have several reproductive systems and life histories. Pest reproductive system can be specified using the repro argument, which accepts "asexual" (haploid), "sexual" (monoecious), and "biparental" (dioecious). For "sexual" pests, the selfing argument specifies if self-fertilisation is (TRUE) or is not (FALSE) allowed. At the start of a simulation, pests are initialised in a random location. Initialised pests are of age zero if rand_age = FALSE or a random age from zero to max_age if rand_age = TRUE. Following initialisation, a single time step proceeds with landscape change (see above), pest aging and metabolism, feeding, pesticide consumption, movement, reproduction, mortality, and immigration (Figure 2). Feeding, pesticide consumption, movement, and reproduction all depend on pest age. Pests feed and consume pesticide from ages min_age_feed to max_age_feed, move from ages min_age_move to max_age_move, and reproduce from ages min_age_reproduce to max_age_reproduce (all inclusive). Food accumulated is lost during aging if baseline_metabolism > 0 and pest age is within min_age_metabolism and max_age_metabolism. The option to set minimum and maximum ages for events makes it possible to model pests with much different life histories (e.g., Sudo et al. 2018).

Figure 2. Overview of simulated events in the resevol R package. Note that metabolism, feeding, pesticide consumption, movement, and reproduction are all subject to a minimum and maximum pest age. Consequently, simulation order might not reflect the order of events from the perspective of a focal pest (e.g., pests might move from ages 1-2, but only feed from ages 2-4). Crops and pesticides are also not necessarily rotated in each time step (see Landscape). Statistics collected within a time step are printed to a CSV file.

In each time step, pests feed in a random order, consuming the crop on their landscape cell. Each pest consumes an amount of crop as specified by the food_consume argument, which takes a vector with as many elements as their are crops (i.e., if crop_number = 2, then food_consume has two elements, the first and second defining consumption of crops 1 and 2, respectively). If crop amount on the landscape cell exceeds pest consumption ability, then pests consume their maximum amount, and this amount is removed from the landscape cell. If crop amount is less than pest consumption ability, then pests consume whatever crop is left and crop amount is reduced to zero. Pesticide consumption works identically to crop consumption, except that the amount of pesticide on a landscape cell is not decreased. Hence, each pest can potentially feed and be affected by the pesticide of their focal landscape cell. Pests simply consume an amount of pesticide as specified by the pesticide_consume argument, which also takes a vector with as many elements as their are pesticides. Pesticide consumption affects pest survival and reproduction using pesticide_tolerated_surv and pesticide_tolerated_repr arguments.

After interacting with their landscape cell, pests can move. Each pest visits a number of cells during movement, as is specified by the parameter movement_bouts. Individual movement bouts occur in a random order across pests. During a movement bout, a pest can travel to any cell within a value defined by movement_distance from their current location, which could include their current location (i.e., moving zero distance). Upon arrival to a cell, a pest can feed if they are of an appropriate feeding age and feed_while_moving = TRUE. The pest also consumes pesticide if they are of the appropriate age and pesticide_while_moving = TRUE. Having pests feed and consume pesticide in a random order while moving among landscape cells can model a population competing for food and encountering pesticides on a shorter time scale than an individual time step.

After all pests finish moving, pests reproduce. Offspring production is possible for asexual, monoecious, or female pests. Pest expected offspring number is defined by a fixed parameter if reproduction_type = "lambda", or is calculated from the amount of food consumed if reproduction_type = "food_based". The former requires specifying lambda_value, which becomes the rate parameter for sampling offspring number from a Poisson distribution. The latter requires specifying a real value for food_needed_repr, which is the amount of food needed to produce one offspring. For food-based reproduction, the total amount of food consumed is divided by food_needed_repr, then floored to determine offspring number. Sexual reproduction requires a mate of reproductive age that is either monoecious or male, and within range of the reproducing focal pest (potential including the focal pest, if selfing = TRUE). A potential mate is within range if it is within an integer number of cells from the focal pest, as defined by mating_distance (e.g., if mating_distance = 0, then mates must share a landscape cell). All available potential mates sire offspring with equal probability, and reproducing pests are assumed to mate multiply (i.e., paternity is a fair raffle for all offspring). If a carrying capacity at birth is set (K_on_birth > 0) and total offspring number in the population exceeds this capacity, then offspring are removed at random until they are within carrying capacity. A real value immigration_rate specifies the rate parameter for Poisson random sampling of the number of immigrants added to the population at the end of a time step. Immigrants are initialised in the same way as pests are at the start of the simulation. Hence, any spatial structure or evolution that occurs during the simulation does not affect immigrant locations, genomes, or traits.

Simulation output

Simulation output is typically large, so output is printed in two CSV files, both of which are created in the working directory. The first file “population_data.csv” prints population level data over time, including population size, mean age, sex ratio, mean food and pesticide consumed of each type, mortality rate, and mean trait values. The second file “individuals.csv” prints all information, including full genomes and traits (columns), for every individual (rows) in the population. The printing of individual level data is disabled by default. It can be turned on for all time steps by setting print_inds = TRUE, but this should be done with caution because it can create extremely large files. Instead, individual level data can be printed for only the final time step by setting print_last = TRUE. Output produced by run_farm_sim is a list of two elements, which includes a vector of parameter values used in the simulation and the final state of the landscape as an array.

Example of individual-based simulations

In resevol v0.3, the run_farm_sim function includes 68 arguments, which specify a wide range of possible simulation conditions affecting landscape and pest characteristics. These arguments are explained the package documentation, and in Supporting Information 2, which demonstrates an advanced case study with a custom landscape and complex pest genomes and life history. Here we focus on a simple simulation with asexually reproducing pests that have three loci and two traits (Figure 1). We use the the genome generated in new_network from Section 2, in which traits 1 and 2 have variances of 1.0008181 and 1.027134, respectively, with a covariance of -0.3925799.

Traits 1 and 2 will define the realised rate of uptake of the two separate pesticides, so we model a system in which there is a potential trade-off for pesticide susceptibility. We use a simple 64 × 64 cell landscape with nine farms. Each farm grows the same crop and uses one of two randomly intialised pesticides, which are rotated every 16 time steps. Hence, we can conceptualise 16 time steps as a single growing season.

In our example, each cell produces four crop units, all of which can potentially be consumed. Pests consume up to one unit of crop on their landscape cell per time step. Pest survival and reproduction is food-based, and pests need to consume one unit of crop by age two to survive and reproduce. Pests initialised at the start of the simulation are randomly assigned an age from 0-4 with equal probability. Pests have a maximum age of four, and they feed in ages 0-2, move up to two cells in ages 3-4, and reproduce in age four. Pests can uptake (i.e., ‘consume’) pesticides in ages 0-2, and if they consume any pesticide, then they will die. In each time step, a mean of 10 immigrants arrive. We simulate 160 time steps (10 growing seasons) using the function below:

sim <- run_farm_sim(mine_output = new_network, repro = "asexual", 
                    pesticide_number = 2, pesticide_init = "random", 
                    pesticide_consume = c("T1", "T2"), farms = 9,
                    pesticide_rotation_time = 16, pesticide_rotation_type = 3,
                    pesticide_tolerated_surv = 0, pesticide_per_cell = 1,
                    crop_rotation_time = 16, crop_number = 1, crop_per_cell = 4, 
                    food_consume = 1, reproduction_type = "food_based", 
                    food_needed_surv = 1, food_needed_repr = 1, max_age = 4, 
                    min_age_feed = 0, max_age_feed = 2, min_age_move = 3, 
                    max_age_move = 4, min_age_reproduce = 4, print_gens = FALSE,
                    max_age_reproduce = 4, age_pesticide_threshold = 2, 
                    rand_age = TRUE, move_distance = 2, immigration_rate = 10, 
                    time_steps = 160, print_last = TRUE, xdim = 64, ydim = 64,
                    trait_means = c(0.1, 0.1), land_edge = "torus");

Any arguments to run_farm_sim not included above are set to default values. Output files can be used to plot ecological and evolutionary dynamics of pests. Figure 3a shows the pest population change both across and within seasons. Due to the specific parameter values chosen, clear patterns reflecting pest cohorts emerge. At the start of a season when the most crop is available, pests eat and abundance increases. As less crop remains in a season, pest abundance decreases. Mean food consumed varies over the course of a season caused by the varying frequency of pests in different life history stages (Figure 3c). In the first 50 time steps, there is some consumption of pesticides 1 and 2 (Figure 3b), but both traits rapidly evolve to negative values making both pesticides ineffective (Figure 3d). In this case, the trade-off in pesticide consumption is not strong enough to maintain susceptibility to either pesticide in the population. Figure 4 shows the location of pests on the landscape in the last time of the simulation (left panel), and how different pesticides are currently being applied on the landscape (right panel). Because the number of crops and pesticides, the number and nature of traits, and the size of pest genomes all vary depending on simulation parameters, no plotting functions are introduced in the resevol package. Instead, methods and code for producing plots such as those in Figures 3 and 4 are explained in Supporting Information 2.

Figure 3. Agricultural pest ecological and evolutionary dynamics over 160 time steps from an individual-based simulation using the resevol R package. Panels show (a) pesticide abundance change, (b) mean realised amount of pesticides 1 and 2 uptaken per pest, (c) mean food consumed per pest, and (d) mean value of evolving traits 1 and 2 underlying pest uptake over time. Note that only pests with positive values for Traits 1 or 2 can uptake Pesticide 1 or 2, respectively. Pests with negative trait values will be unaffected by corresponding pesticides (i.e., pesticide consumption is a threshold trait), hence the difference between realised pesticide consumption (b) and the traits underlying it (d). White and grey vertical stripes indicate seasons of a single crop and pesticide application.

Figure 4. Locations of pests (black) across a landscape that includes nine farms (coloured blocks) in the last time step of a simulation using the resevol R package (left panel). The right panel shows which farms apply pesticide 1 (dark grey) and 2 (light grey).

Conclusions

Insecticide resistance to pesticides is a widespread problem that affects global food security (Tabashnik et al. 2014; Bradshaw et al. 2016; Lykogianni et al. 2021). Many models have been developed to investigate the evolution of resistance under different ecological and evolutionary conditions (e.g., Curtis 1985; Mani 1985; Ives and Andow 2002; Liang et al. 2013; Liang, Tang, and Cheke 2018; Haridas and Tenhumberg 2018; South and Hastings 2018; Sudo et al. 2018; Saikai, Hurley, and Mitchell 2021). The resevol R package makes it possible to rapidly develop and simulate myriad individual-based models of resistance evolution. It also introduces a novel approach to modelling complex pest genetic architecture, using an evolutionary algorithm to generate haploid or diploid loci that map to pest traits with pre-specified covariances (Figure 1; Supporting Information 1). Agricultural landscapes and pest life histories are highly customisable, allowing targetted models that can simulate specific real-world case studies (Supporting Information 2). The breadth of possible models that can be simulated with resevol also makes it a useful tool for developing theory on pest management, and even more generally on the evolution and ecology of individuals with complex traits on a heterogeneous landscape.

Acknowledgements

This software was developed as part of the project for Enhancing Diversity to Overcome Resistance Evolution (ENDORSE). The ENDORSE project is a joint Newton funded international partnership between the Biotechnology and Biological Sciences Research Council (BBSRC) in the UK and the São Paulo Research Foundation (FAPESP) in Brazil under BBSRC award reference BB/S018956/1 and FAPESP award reference 2018/21089-3. ENDORSE is a partnership among Universidade Estadual Paulista (UNESP), the University of Stirling (UoS), and the Centre for Agricultural and Biosciences International (CABI).

Authors contributions

All authors contributed to planning software development and writing. ABD wrote and maintains package code. CRM contributed code for data visualisation.

Availability

The resevol R package can be downloaded from CRAN (https://cran.r-project.org/package=resevol) or GitHub (https://bradduthie.github.io/resevol/). The package is open source under GNU Public License.