Individual-based modelling of the efficacy of fumigation tactics to control lesser grain borer (Rhyzopertha dominica) in stored grain

Journal of Stored Products Research 51 (2012) 23e32

Contents lists available at SciVerse ScienceDirect

Journal of Stored Products Research journal homepage: www.elsevier.com/locate/jspr

Individual-based modelling of the efficacy of fumigation tactics to control lesser grain borer (Rhyzopertha dominica) in stored grain Mingren Shi a, c, *, Patrick J. Collins c, e, James Ridsdill-Smith b, c, d, Michael Renton a, c, d a

M084, School of Plant Biology, FNAS, The University of Western Australia, 35 Stirling Highway, Crawley, WA 6009, Australia M092, School of Animal Biology, FNAS, The University of Western Australia, 35 Stirling Highway, Crawley, WA 6009, Australia c Cooperative Research Centre for National Plant Biosecurity, Australia d CSIRO Ecosystem Sciences, Underwood Avenue, Floreat, WA 6014, Australia e Agri-Science Queensland, Department of Agriculture, Fisheries and Forestry, Ecosciences Precinct, GPO Box 267, Brisbane, QLD 4001, Australia b

a r t i c l e i n f o

a b s t r a c t

Article history: Accepted 8 June 2012

Increasing resistance to phosphine (PH3) in insect pests, including lesser grain borer (Rhyzopertha dominica) has become a critical issue, and development of effective and sustainable strategies to manage resistance is crucial. In practice, the same grain store may be fumigated multiple times, but usually for the same exposure period and concentration. Simulating a single fumigation allows us to look more closely at the effects of this standard treatment. We used an individual-based, two-locus model to investigate three key questions about the use of phosphine fumigant in relation to the development of PH3 resistance. First, which is more effective for insect control; long exposure time with a low concentration or short exposure period with a high concentration? Our results showed that extending exposure duration is a much more efficient control tactic than increasing the phosphine concentration. Second, how long should the fumigation period be extended to deal with higher frequencies of resistant insects in the grain? Our results indicated that if the original frequency of resistant insects is increased n times, then the fumigation needs to be extended, at most, n days to achieve the same level of insect control. The third question is how does the presence of varying numbers of insects inside grain storages impact the effectiveness of phosphine fumigation? We found that, for a given fumigation, as the initial population number was increased, the final survival of resistant insects increased proportionally. To control initial populations of insects that were n times larger, it was necessary to increase the fumigation time by about n days. Our results indicate that, in a 2gene mediated resistance where dilution of resistance gene frequencies through immigration of susceptibles has greater effect, extending fumigation times to reduce survival of homozygous resistant insects will have a significant impact on delaying the development of resistance. Ó 2012 Elsevier Ltd. All rights reserved.

Keywords: Individual-based model Two-locus simulation Phosphine resistance Lesser grain borer Management tactics

1. Introduction The safe storage and supply of cereal grains and foods depends on control of potentially highly destructive insect pests, particularly in tropical and subtropical regions. Owing to its international market acceptance and a lack of acceptable, cost-effective alternatives, disinfestation with phosphine (PH3) fumigant is a fundamental tool used world-wide in the management of these insects. However, the development of resistance to phosphine by the lesser grain borer, Rhyzopertha dominica, a serious cosmopolitan pest of * Corresponding author. M084, School of Plant Biology, FNAS, The University of Western Australia, 35 Stirling Highway, Crawley, WA 6009, Australia. Tel.: þ61 8 6488 1992; fax: þ61 8 6488 1108. E-mail address: [email protected] (M. Shi). 0022-474X/$ e see front matter Ó 2012 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.jspr.2012.06.003

stored cereal grains, seriously threatens effective insect pest management (Collins, 2006; Emekci, 2010). The response of pest managers to resistance has been to increase phosphine concentrations and exposure periods (Collins et al., 2005) and, in cases where control cannot be achieved, apply residual insecticides. The latter is a last resort, however, as it can limit market access. Other measures adopted to combat resistant insects include strategic use of sulfuryl fluoride, intensive storage hygiene and increased insect population monitoring and resistance testing (Nayak et al., 2010). A strategy relying primarily on a single fumigant backed by a limited number of alternatives is highly risky, however, and will require very careful management to sustain. The development of resistance to pesticides in insects is affected by a variety of interacting influences, including genetic, biological/ ecological and management (operational) influences (National

24

M. Shi et al. / Journal of Stored Products Research 51 (2012) 23e32

Research Council, 1986). Many factors can affect whether an insect population will survive under phosphine fumigation treatments. Some of the more important factors include gas concentration, duration of fumigation and developmental stages of the insects present. Our research is aimed at contributing to the development of viable, long-term strategies to support the management of phosphine resistance. Computer simulation models can provide a relatively fast, safe and inexpensive way to project the consequences of different assumptions about resistance and to weigh the merits of various management options. We used stochastic individual-based modelling that explicitly represents the fact that R. dominica populations consist of individual beetles, each of a particular genotype and a particular life stage. Individual-based approaches are relatively easily adapted to incorporate the biological attributes that we want to investigate, such as different initial frequencies of genotypes and different proportions of life stages. Resistance to phosphine is an inherited trait and two major genes are responsible for the strong phosphine resistance in R. dominica (Collins et al., 2002). These two genes act in synergy to cause a significantly increased resistance to phosphine compared with either one of the resistance genes on their own (Schlipalius et al., 2008). Our previous research (Shi et al., in press) showed the importance of basing resistance evolution models on realistic genetics and that using simplified one-locus models to develop pest control strategies runs the risk of not correctly identifying tactics to minimise the incidence of pest infestation. Hence, our simulations were carried out using the two-locus model. In this paper we used our individual-based two-locus model to address three questions about the management tactics of single phosphine fumigation by investigating some biological and operational factors that influence the development of phosphine resistance in R. dominica. First, which application tactic is more effective: long exposure time with low concentration, or short exposure time with high concentration? This involved looking at the impact on insect life stage (egg, larva, pupa, adult), each of which has an inherently different tolerance to phosphine, on the efficacy of two equi-toxic phosphine treatments: Fum1: “0.0685 mg/l  16 days (exposure for 16 days at a concentration of 0.0685 mg/l)” and Fum2: “0.2 mg/l  8 days (exposure for 8 days at a concentration of 0.2 mg/ l)”. The second question is: within an insect population of a fixed total size, what is the impact of different initial genotype frequencies on the efficacy of phosphine fumigation? Following on from this, how long should fumigation time be extended to provide an equally effective level of control, if the initial frequency of strongly resistant beetles is n times greater than the original one? This could represent populations with different histories of exposure to phosphine fumigation. Third, given a fixed initial genotype frequency, how does the level of infestation, i.e. the total number of insects within a grain storage, impact on the effectiveness of phosphine treatment? This may represent varying levels of grain hygiene maintenance or varying amounts of movement of pests from infested places outside to inside a storage facility.

2.1. Assumptions regarding genotypes and resistance Resistance to phosphine is an inherited trait. Our simulation model was constructed based on results from Collins et al. (2002) and Schlipalius et al. (2002). Their research revealed that the combination of alleles at two loci, rph1 and rph2, confers strong resistance while rph1 by itself is responsible for the weak resistance phenotype. It seems that both rph1 and rph2 individually express a relatively low level of resistance but when they occur in the same insect the resistance mechanisms synergise, producing a much higher level of resistance. Based on this data, we assume there are two possible alleles on each of these two loci, a susceptibility allele and an incompletely recessive resistance allele. At both loci, the susceptibility allele is assumed to be relatively common initially, and the resistance allele to be relatively rare. At each locus there are thus three possibilities, which we denote as s (homozygous susceptible), h (heterozygous) and r (homozygous resistant). Our two-locus simulation model thus includes nine possible genotypes, which can be denoted as ss, sh, sr, hs, hh, hr, rs, rh and rr. For example, ss denotes the genotype with both loci homozygous susceptible, and rh denotes the genotype with the first locus homozygous resistant, and the second locus heterozygous. 2.2. Overview of individual-based model dynamics and assumptions The simulated dynamics for each individual beetle at each daily time step during the simulation are illustrated in Fig. 1, and the default parameter values used in the model are listed in Table 1 with a brief description. The life stages include egg, larva, pupa and adult. We separate the adult stage into two: adult 1 (immature unable to lay eggs) and adult 2 (mature able to lay eggs) in simulations but we merge the counts of adult 1 and adult 2 into a single adult stage in the results. Life history parameters of each life stage were estimated from published experimental data (Collins et al., 2002, 2005; Daglish, 2004) based on the assumed temperature and relative humidity. A number of processes within the simulation are determined stochastically:  The ‘time remaining within current life stage’ (TRICLS) for each individual beetle is drawn randomly from a normal distribution, with mean and standard deviation depending on the life stage (Table 1) in “Initialize” or “Enter the next stage” step (Fig. 1).  Whether the beetle survives through the day during a period of fumigation is drawn from a Bernoulli distribution with survival probability that depends on the genotype of the beetle, and the concentration and exposure time of the fumigant, as explained below.  The sex of the egg is drawn from a Bernoulli distribution, with an even probability of being male or female (male eggs removed from simulation). Initialise: genotype

2. Model and methods An overview of the two-locus model is given here; full details of the model have been provided previously (Shi and Renton, 2011; Shi et al., in press, Shi et al. 2012). The model assumptions regarding genotypes and resistance are described first, followed by explanation of overall model dynamics and simulation processes. Then we describe the methods we used to estimate the two kinds of finite daily survival rates under the two PH3 treatments. Finally, we describe in detail how we investigated the three central questions posed in the Introduction.

life stage TRICLS

To the next day TRICLS

1

Determine: No

TRICLS <= 0? life [ Current stage ends? ] Yes

Enter the next stage Set new TRICLS

Yes

Dead

No

Survival? Yes

Yes

Is life time over?

number of eggs genotype / sex

No

currently No under fumigation?

egg duration Yes

Able to lay eggs?

No

Fig. 1. The overall model dynamics occurring each day for each individual beetle during a simulation.

M. Shi et al. / Journal of Stored Products Research 51 (2012) 23e32

leakage (Banks, 1989) is negligible. Temperature and relative humidity within the store remains the same so that life stage durations are constant.

Table 1 Default parameters, variables and abbreviations for our two-locus model. Parameter Description

Value megg ¼ 11.9, mlarva ¼ 36.5, mpupa ¼ 9.6, madult1 ¼ 15, madult2 ¼ 102 Segg ¼ 1.5, Slarva ¼ 4.6, Spupa ¼ 1.2, Sadult1 ¼ 0, Sadult2 ¼ 15 0.2242 100,000 (100 K) Variable Variable

mi

Mean number of days at life stage i

si

Standard deviation of number of days at life stage i

bD N0 TRICLS TSR

Daily (finite) birth rate Starting number of (female) beetles Time remaining in current life stage Total survival rate at the end of a fumigation treatment Cumulative survival rate each day Variable during a fumigation period Geometrical average daily survival rate Variable Different daily survival rate Variable

CSR GADSR DDSR

25

2.3. Estimation of survival rates under fumigation We used probit models, as we found they had smaller least squares (numerical) errors for our data than logistic and Cauchy models (Shi and Renton, unpublished results), to predict the total mortality or survival (¼ 1 e mortality) rates of adult beetles under different concentrations and exposure times. We used the fourparameter probit model

Y ¼ a þ b1 logðtÞ þ b2 logðCÞ þ b3 logðtÞlogðCÞ

(1)

for the beetles of genotypes hs, rs and rr, as it provided a better fit of the available data than the 2- and 3-parameter models in terms of smaller least squares (numerical) errors. Here, Y is the probit (the inverse cumulative distribution function) value of mortality, C is concentration (mg/l) and t is exposure time (hours). We used the two-parameter probit model

 The number of eggs produced by an individual in a day is determined by drawing randomly from a Poisson distribution with mean equal to the daily birth rate bD (Table 1).  The genotype of each new egg is determined by drawing randomly from a multinomial distribution based on the maternal and paternal genotype and the offspring genotype table (Shi and Renton, 2011).

Y ¼ a þ b logðCtÞ

(2)

to predict the mortality rates for other genotypes, because t is a constant (48 h) in the relevant data set. These models were fitted to data on final mortality rates for different genotypes, concentrations and durations (Collins et al., 2002, 2005; Daglish, 2004) and the fitted parameters for the two models are listed in Table 2. Directly substituting the fixed C value and the series of t values 1, 2, ., T into Eqs (1) and (2) results in the predicted cumulative survival rates (CSR) after 1, 2, ., T days. The CSR for the final Tth day is the total survival rate (TSR) for the full C  T treatment. We wanted to design two phosphine treatments, Fum1 and Fum2, to test whether long exposure time with low concentration, Fum1, or short exposure time with high concentration, Fum2, is more effective, i.e. results in lower survival rates, taking into account that R. dominica tolerance to phosphine varies with life stage. We wanted these two treatments to be equi-toxic, meaning that for a population of only adult beetles they would result in the same total mortality. We chose Fum2 to be the treatment “0.2 mg/ l  8 days” (C ¼ 0.2 mg/l and t ¼ 8 days). Under this Fum2 treatment, the total survival rate of the rr and rh beetles are 0.007888 and 7.2  1010 respectively, and <1.04  1021 for all other genotypes. In other words, this treatment kills about 99.2% (z1  0.007888) of the rr beetles and effectively eradicates the other genotypes. This ‘threshold’ treatment achieves a high kill rate but may allow a small fraction of highly resistant beetles to survive. Achieving this dose for the whole duration throughout the storage facility can be seen as an ideal target fumigation strategy. In reality doses achieved in storage facilities are likely to be more variable in space and time (also see Discussion).

Full details of model dynamics and parameter values are included in Shi et al. (in press). A 1:1 sex ratio was assumed, which allowed us to simulate the female beetles only, with the assumption that total number and the genotype frequency for the male beetles was the same as for the females. The rate of development and survival of various life stages of the lesser grain borer depend on several environmental factors, among them temperature and moisture content of the grain (or equilibrium relative humidity) were the main factors considered in our simulations. For this model, we assumed grain type wheat stored at a typical field temperature (T) of 25  C (Cassells et al., 2003) and relative humidity (r.h.) of 70%. The latter was chosen as much of the published life history data are provided at 70% r.h. (Shi et al. 2012). In addition, these environmental conditions are recommended by FAO for bioassay (FAO, 1975) and used by many researchers in their experiments (e.g. Collins et al., 2000, 2002, 2005; Daglish, 2004; Herron, 1990; Pimentel et al., 2007; Schlipalius et al., 2008). It is assumed that no eggs are laid during periods of fumigation. Mating occurs randomly; for example, resistance genetics does not affect choice of mate. The grain food supply is non-limiting and its quality does not affect the natural birth rate of the beetle. Fumigant concentrations do not vary with time or location within the storage facility or silo. Variation with time or space due to uptake or release of gas from or into grain (sorptionedesorption), diffusion, or Table 2 The fitted parameters of two- and four-parameter probit models. Genotype

ss hs rs sh hh rh sr hr rr

Two-parameter probit model

Four-parameter probit model

a

b

a

b1

b2

b3

15.032386

9.229083 11.284676 10. 413046

3.776399 15.575413

6.964954 0.047656

1.010451 4.701759

12.232356

10.386287

3.101974

1.190773

10.854928 7.954711 1.682448 7.043248 3.944565

5.913329 5.913329 5.913329 5.913329 5.913329

26

M. Shi et al. / Journal of Stored Products Research 51 (2012) 23e32

To set a treatment with a lower concentration C and longer exposure time t to achieve the same 99.2% total kill rate of the adult rr beetles we doubled the exposure time, i.e. setting t ¼ 16 days, and then the C value required to make the treatment equi-toxic could be estimated by rearranging Eq. (1):

logðCÞ ¼ ½Y  a  b1 logðtÞ=½b2 þ b3 logðtÞ

(3)

with t ¼ 16  24 (hours), Y ¼ 7.414 (corresponding to mortality 1  0.007888), and the parameters listed in Table 2 for the rr genotype. The result is C ¼ 0.0685 mg/l. Thus we had two equi-toxic fumigation treatments:

Fumigation 1ðFum1Þ : “0:0685 mg=1  16 days”; Fumigation 2ðFum2Þ : “0:2 mg=1  8 days”:

(4)

As desired, these treatments predicted the same TSR for the adult rr beetles. It can be seen from Table 3 that the CSR after 8 days under Fum2 is indeed exactly equal to the CSR after 16 days under Fum1, and in fact the nth day’s CSR under Fum2 is very close to the (2n)th day’s CSR under Fum1 in general. The probit models predicted total survival rate (TSR) from a fumigation treatment lasting a number of days, but our individual-based simulation needed predictions of the daily mortalities within the fumigation period (the finite daily survival rate). There are two ways to estimate the finite daily survival rate under a PH3 treatment C  T at a fixed concentration C (mg/l) and a range of times (1, 2, ., T days). The first way to estimate the finite daily survival rate is based on the assumption that the survival rate is the same each day during a fumigation treatment, in which case the finite daily survival rate of genotype x is [TSR(x)]1/T. We can call this geometrical average daily survival rate (GADSR). The estimated GADSRs for the nine genotypes under the two treatments Fum1 and Fum2 are listed in Tables 4 and 5. Note that if we use GADSR to design Fum1, just simply set GADSRFum1 ¼ (GADSRFum2)1/2. In other words, for the rr beetles the CSR for one day under Fum2 is exactly equal to two-day’s CSR under Fum1 or the square of the GADSR under Fum1. The second way is to estimate a different daily survival rate (DDSR) for each day: convert the ith day’s CSRi and the previous (i  1)th day’s CSR i  1 into the ith day’s DDSRi(x) (for genotype x) by setting DDSR1(x) ¼ CSR1(x) and then letting DDSRi(x) ¼ CSRi(x)/ CSRi  1(x), for i ¼ 2, ., T. (the current day’s DDSR is equal to the ratio of the current day’s CSR to the previous day’s CSR). The DDSR of the rr beetles under Fum1 and Fum2 thus derived are listed in Table 6, and illustrated in Fig. 2. If the simulation is concerned only with the results after fumigation, then we can use the former approach with equal daily survival rates (GADSR), as it is simpler. If the simulation is also concerned with accurately representing what happens during the fumigation process then we should use the latter approach (DDSR), even though it is more complex and will increase the simulation Table 3 Cumulative survival rates (CSR) under Fum1 up to the (2n)th day and Fum2 up to the nth day for the adult rr beetles (ratio ¼ [Fum2 at (n)th day]/[ Fum1 at (2n)th day] for rr beetle at any stage). n

1

2

3

4

Fum1/(2n)th day Fum2/(n)th day Ratio

0.99999999 1. 1.00000001

0.99867 0.99958 1.00091

0.922132 0.951103 1.031417

0.6160571 0.6779515 1.1004686

n

5

6

7

8

Fum1/(2n)th day Fum2/(n)th day Ratio

0.28196 0.32136 1.13973

0.09859 0.11116 1.12749

0.02924 0.03144 1.07548

0.007888 0.007888 1.000000

The bold values represents the total survival rates under Fum1 (16 days) and under Fum2 (8 days) are the same 0.007888.

Table 4 The geometrical average daily survival rate (GADSR) for the nine genotypes for Fum2 treatment (0.2 mg/l  8 days). 1st gene

s h r

2nd gene s

h

r

<1.  1011 <1.  1011 0.000304

<1.  1011 <1.  1011 0.007196

6.33  1010 0.00238 0. 54591

time. We used DDSR for the simulation of the impact of different life stages on the effect of fumigation since we wanted to see the details of what happens during the fumigation process in each life stage. But we used GADSR for the other two simulations as we were concerned only with the survival numbers after fumigation. Full details for the estimation of survival rates are included in Shi and Renton (2011) and Shi et al. (2012). 2.4. Question 1: which is more important: concentration or exposure time? Flinn et al. (2001) reported that the days spent to reach 95% mortality exposed to 180 ppm (0.25 mg/l) at 25  C for life stages of R. dominica were as follows:

Egg; ð1:0;

Larva; 0:4;

Pupa; 0:5;

Adult 0:3Þ

(5)

For simplicity, we assume that the relative tolerance of stages does not change with genotype. We also assumed that the relative tolerance to PH3 of each life stage and therefore the TSR (total survival rate) was equal to the ratio of the times to 95% mortality provided in Eq. (5). The predicted time to 95% mortality of adults exposed to 0.0685 mg/l is approximately 12.5 days. If we regard this as the TSR for adults then the DDSR of the rr beetles at other life stages under Fum1 for each day are:

DDSRE ¼ ð1:0=0:3Þ1=12:5 DDSRA z1:1055 DDSRA ; DDSRL ¼ ð0:4=0:3Þ1=12:5 DDSRA z1:0243 DDSRA and DDSRP ¼ ð0:5=0:3Þ1=12:5 DDSRA z1:0435 DDSRA

(6)

where DDSRA is the daily survival rate of the rr beetles at adult stage, and DDSRE DDSRL and DDSRP are the daily survival rate of the rr beetles at the egg, larva and pupa stage, respectively. We can now set the DDSR of the rr beetles at other life stages under Fum2 so that the TSRs under Fum1 and Fum2 at each stage are almost the same. That is, under Fum2,

DDSRE z1:2222 DDSRA ; DDSRL z1:0491 DDSRA and DDSRP z1:0889 DDSRA

(7)

We emphasize again that the beetles of other genotypes are all killed under Fum1 and Fum2. We ran the model starting with equal proportions of the four life stages PS ¼ (0.25, 0.25, 0.25, 0.25). We also started with a resistance allele frequency of q ¼ v ¼ 0.62947, where q and v are the frequency of the resistance alleles at the 1st and 2nd locus respectively, corresponding to a high rr genotype frequency of 0.157, representing Table 5 The geometrical average daily survival rate (GADSR) for the nine genotypes for Fum1 treatment (0.0685 mg/l  16 days). 1st gene

s h r

2nd gene s

h

r

< 1.  1011 6.3  106 0.02223

< 1.  1011 3.04.  105 0.3801

0.00074 0.12849 0.7389

M. Shi et al. / Journal of Stored Products Research 51 (2012) 23e32

27

Table 6 Different daily survival rates for the nth day (n/DDSRn) of the adult rr beetles under Fum1 (16 days) and Fum2 (8 days). Fum1 Fum2

1/1. 9/0.705 1/1.

2/0.99999999 10/0.649 2/0.9996

3/0.99998 11/0.607 3/0.952

a population on the threshold of exhibiting a serious resistance problem (see Shi and Renton, 2011 for details on setting initial genotype frequencies). We started to fumigate immediately and ran the model for 140 days (20 weeks) in total, enough time for populations to build up again after the fumigation. 2.5. Question 2: impact of initial genotype frequencies on fumigation efficacy To simulate the impact of different initial frequencies of genotypes on the effectiveness of fumigation, we first set up nine sets of initial frequencies of genotypes. We defined the resistance allele frequencies for the nth set (n ¼ 1, 2, ., 9) as qn ¼ vn ¼ (0.1 n)1/4. We then calculated the equilibrium genotype frequencies for the nine sets, resulting in initial frequencies for the rr beetles set of:

0:1; 0:2; 0:3; .; 0:8; 0:9 That is, the initial frequency of the nth set for rr beetles is n times the original one (0.1). We set the initial proportions of the four life stages to be equal PS ¼ (0.25, 0.25, 0.25, 0.25). We started to fumigate immediately under Fum2 and ran the model for 140 days (20 weeks) in total. We also ran the model to fumigate from the 1st day to 40th day at the Fum1 concentration 0.0685 mg/l (and to fumigate from the 1st day to 30th day at the Fum2 concentration 0.2 mg/l to learn how many extra days are required to reach the same control level). Note that to separate the influence of the different factors we ignored the impact of different life stages in this simulation, that is, the GADSRs listed in Table 4 or Table 5 were used for all life stages. The start population number of (female) beetles for this simulation was also 100 K.

4/0.9986 12/0.576 4/0.713

5/0.984 13/0.553 5/0.474

6/0.938 14/0.536 6/0.346

7/0.860 15/0.524 7/0.283

8/0.777 16/0.515 8/0.251

2.6. Question 3: impact of initial population size on fumigation efficacy This simulation was undertaken to test the impact of a range of sizes of initial insect populations within the grain storage, representing various levels of storage hygiene (amounts of grain and other residues around the storage facility that harbour insects and provide sources of initial infestation). Initial genotype frequencies were the same as for Question 1. The initial proportions of the four life stages and the durations of fumigation and simulation were the same as for Question 2. We simulated the Fum1 treatment only, and considered five different initial population numbers: 100 K, 200 K, 400 K, 800 K and 1,600 K. 3. Results 3.1. Question 1: which is more important: concentration or exposure time? The numbers of the rr beetles in all life stages decline at different rates under the two fumigation treatments, and the proportions in each life stage also vary in different ways (Fig. 3). These differences then affect how the populations recover after fumigation ends (Fig. 4). The number of rr beetles after fumigation under Fum2 are about 1.4 times as many as those under Fum1 (Table 7), and this difference carries through to a similar difference at the end of the whole simulation period. The numbers and proportions under Fum1 for the first eight days and those under Fum2 for the first four days at each stage are very similar (left parts of subplots in Fig. 3). In this initial period, the numbers of rr beetle at all life stages decrease. The adult beetles, the least tolerant life stage, decrease most quickly, but there is little difference between the life stages. After this initial period, the differences between the four life stages and the two treatments become more evident (right part of Fig. 3). The decline in numbers under Fum1 is faster than under Fum2, with pupa and egg numbers dropping to zero in Fum1 but not Fum2. 3.2. Question 2: impact of initial genotype frequencies on fumigation efficacy

Fig. 2. The different daily survival rates (DDSR, the curve on the top) and its corresponding cumulative survival rate (CSR from DDSR, having the same start point as the DDSR curve; the 1st day’s DDSR is equal to the 1st day’s CSR), and the geometrically average daily survival rate (GADSR, horizontal line) and the corresponding cumulative survival rate (CSR from GADSR, also having the same start point as the GADSR line) of the adult rr beetles under Fum2 (0.2 mg/l  8 days). Note that the DDSR changes daily, while the GADSR is constant, but both result in the same total survival/mortality (their corresponding CSR curves end at the same point).

As the initial frequency of the rr beetles is increased, the numbers of rr beetles at the end of fumigation and after simulation both increase as well, as shown by the numbers in Table 8 and the fact that the two curves in Fig. 5 (a) are both increasing. Furthermore, the two curves in Fig. 5(a) are nearly parallel, which indicates that the number of beetles at the end of the simulation is almost a constant multiple of the number of beetles following fumigation. The nine curves in Fig. 5(b) are also parallel, indicating that the effect of different initial genotype frequencies is a simple constant multiplicative one. Each of the nine curves in Fig. 5(b) has two parts; the first part, which is decreasing, is during the fumigation and the second part, which is increasing, is after the fumigation. The relationship between the initial frequency of the rr genotype and the numbers of beetles at the end of fumigation and after simulation is approximately linear (Table 8, Fig. 5(a)), and this is actually true at any time during the simulation (Fig. 5(b)). It can be seen from Table 9 that if the initial frequency is increased to n times

28

M. Shi et al. / Journal of Stored Products Research 51 (2012) 23e32

A

B

C

D

Fig. 3. Numbers (a, b) and proportions (c, d) of beetles with each life stage for the fumigation period under the two treatments Fum1 (a, c) and Fum2 (b, d). The vertical line x ¼ 8 for Fum1 and the line x ¼ 4 for Fum2 separate the figures into two parts, emphasising the fact that the four curves (for four stages) on the left parts of (a) and (b) or of (c) and (d) have very similar behaviour, but the right parts differ. (Note that when numbers reach zero they are not shown because of the log scale).

the original, then we need to extend fumigation duration approximately n days under the Fum1 concentration to achieve a similar control level, whereas under the Fum2 concentration we need to extend fumigation duration approximately a half of n days. 3.3. Question 3: impact of initial population number on fumigation efficacy It can be seen by comparing Fig. 5 and Fig. 6 that the results of this simulation are very similar to those obtained from the previous simulation for Question 2; the two curves in Fig. 6(a) are both increasing and nearly parallel, the five curves in Fig. 6(b) are also parallel, and the ratio of the numbers at the end of fumigation and at the completion of simulation are very similar; if the starting population number is increased n times, the numbers at the end of fumigation and at the completion of simulation also increase about n times (Table 10). Note that the original frequency of rr beetles for this question is 0.157 while that for Question 2 is 0.1; this is why the numbers of rr beetles at the end of fumigation and at the completion of simulation for the two questions differ, even though the starting population number of both was 100 K. Also note that, as our model is stochastic, for any particular model run the increase will only be approximately n times, but when averaged over a large number of model runs, the ratio will be exactly n times.

4. Discussion In this research we used our individual-based, two locus resistance model to investigate three key questions about the impact of a single application of phosphine on control of phosphine-resistant R. dominica populations. As the fumigation parameters of exposure time and concentration are usually under operational control, we asked: at equi-toxic dosages, which is better for insect control; a longer fumigation at lower concentration or a shorter fumigation at a higher concentration? Our model showed that the former is more effective, that is, will kill a higher proportion of insects than the latter. The reason for this is that under a longer exposure period, the more tolerant immature stages (eggs and pupae) have sufficient time to develop to less tolerant stages (larvae and adults) (Fig. 3). As there is no evidence available to indicate that exposure to phosphine may influence immature development time in R. dominica, we have assumed in our model that this does not occur. Phosphine has been shown to delay egg hatch in some insect pests of stored products (Rajendran, 2000; Nayak et al., 2003) but not others (Price and Bell, 1981; Pike, 1994). It is already known that extending the fumigation period (while lowering concentration) will increase toxicity of phosphine (Winks, 1985; Chaudhry, 2000; Collins et al., 2005; Mills and Athie, 1999; Bond, 1984). High concentrations may not increase toxicity and, in reality, they may cause insects to go into a protective ‘narcosis’. It

M. Shi et al. / Journal of Stored Products Research 51 (2012) 23e32

A

B

C

D

29

Fig. 4. Numbers (a, b) and proportions (c, d) of beetles with each life stage for the whole simulation period (fumigation and recovery) under the two treatments Fum1 (a, c) (Note that when numbers reach zero they are not shown because of the log scale.).

appears that the toxic effects of phosphine accumulate slowly in R insects and that the resistance mechanism can be overwhelmed during long exposure periods. These effects are not accounted for in this study, and provide additional support for a long duration strategy, over and beyond the effects due to life stage transitions that we found here. The second key question we asked was: what impact does initial resistance genotype frequency have on phosphine efficacy? In addition, having shown that extending fumigation time is superior to increasing concentration, we asked how long fumigation time should be extended to provide greater control at higher resistance gene frequencies. Our results show that the rate of survival of the pest population will increase proportionally with initial resistance genotype frequency, that is, by about n times if the initial frequency of the rr genotype is increased n times. This is because we assumed the same survival and development rates for all nine possible genotypes, including rr. Thus, to achieve a similar level of control of a population with an initial frequency of the rr genotype that is n Table 7 Total survival numbers (TSN) of the rr beetles after fumigation and at the end of the whole simulation (fumigation and recovery) under the Fum1 and Fum2 scenarios. Treatment

Fum1 (0.0685 mg/l  16 days) Fum2 (0.2 mg/l)  8 days)

Total survival number (TSN)

times higher than another population with a low initial frequency, we need only increase the fumigation time by approximately n days rather than having to multiply the original fumigation duration by n times. This again shows that increased duration of fumigation is an efficient way to increase the efficacy of fumigation. Our third key question was: what impact does the level of infestation or total number of insects within a grain storage have on the effectiveness of a phosphine fumigation? Obviously, the number of rr insects surviving a fumigation will increase about n times if its initial population size is increased by n times. In addition, the initial number of each genotype is always the product of two factors: the total initial population number and the initial frequency of the genotype. If we fix one of these two factors and increase the other factor by n times the original one, this always

Table 8 The numbers of the rr beetles at the end of Fum1 fumigation (EF) and at the end (or completion) of the whole simulation (ES) for each of the initial genotype frequencies of the rr beetles (f(rr)) tested. For each genotype frequency, the ratio of the EF and ES results relative to the EF and ES results at the lowest initial frequency of 0.1 are also given. f(rr)

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

89 1.

139 1.56

225 2.53

324 3.64

395 4.44

444 4.99

507 5.70

620 6.97

723 8.12

1417 1.

2173 1.53

3599 2.54

5328 3.76

6271 4.43

6834 4.82

8132 5.74

9577 6.76

11,593 8.18

After fumigation

After simulation

EF Ratio

160 224

2259 3077

ES Ratio

30

M. Shi et al. / Journal of Stored Products Research 51 (2012) 23e32

A

A

B

B

Fig. 5. The number of rr beetles at the end of the fumigation period and at the completion of running the whole simulation (fumigation and recovery) for each of the 9 initial genotype frequencies (a), and the number of rr beetles over time during the whole simulation (b).

results in an initial number of each genotype equal to n times the original starting number of that genotype. To answer the second key question discussed previously, we fixed the starting population number at 100 K but varied the initial frequency of the rr genotype beetles. To answer this third question, we fixed the initial frequency but varied the initial population number. The results are simplified by the fact that our two PH3 treatments killed all beetles of all genotypes besides rr. Hence the simulations we used to answer key questions 2 and 3 are quite alike, and gave very similar results during the fumigation period (Figs. 5 and 6). Moreover, the daily birth rates in the two simulations were the same, resulting in very similar behaviours after fumigation. Thus we can draw a very similar conclusion for the third question as for the second: to achieve a similar level of control of a population that is n times Table 9 The number of extra days exposure needed under Fum1 and Fum2 to achieve a similar control level to the first simulation when the rr genotype frequency is 0.1, for various cases where the initial rr genotype frequency is higher that 0.1. Initial frequency of rr Extra days under Fum1 Extra days under Fum2

0.2 3 1

0.3 4 2

0.4 6 3

0.5 6 3

0.6 7 3

0.7 7 4

0.8 8 4

0.9 8 4

Fig. 6. The numbers of rr genotype beetles at the end of the fumigation period and at the end of the whole simulation (a), and over time during the whole simulation period (b), for initial populations of 100 K, 200 K, 400 K, 800 K and 1600 K insects, given a fixed initial genotype frequency and equal initial proportions of the four life stages.

larger than another population, we need only to increase the fumigation time by approximately n days rather than having to multiply the original fumigation duration by n times. This once again shows that increased duration of fumigation is an efficient way to increase the efficacy of fumigation. Our two- and four-parameter probit models (1) and (2) fit the experimental data sets of Collins et al. (2002, 2005) and Daglish (2004) very well and this allowed us to accurately predict mortality of R. dominica under a variety of phosphine treatments with different concentrations and exposure times (Shi and Renton, 2011). In previously published resistance modelling research “survivorship was not explicitly included in the model because adequate data were not available” (Hagstrum and Flinn, 1990), and when different mortalities for different genotypes were included, they were only roughly divided into a few levels (e.g. Tabashnik, 1989; Longstaff, 1988) and a simplified single gene model was used. In other cases, mortalities were varied with temperature and moisture in some detail but differences due to concentration, exposure time, or genotype were not included (e.g. Flinn et al., 1992). To our knowledge, no previous models have included mortality predictions that vary with concentration, exposure time,

M. Shi et al. / Journal of Stored Products Research 51 (2012) 23e32 Table 10 The number of the rr beetles at the end of the fumigation period and at the completion of the whole simulation for the five initial population numbers, and the corresponding ratios of numbers resulting from higher starting population numbers to the numbers resulting from a starting population of 100 K insects, representing the lowest starting number. Start population number

100 K

200 K

400 K

800 K

1600 K

At the end of fumigation Ratio

126 1.

235 1.87

459 3.64

1045 8.29

1992 15.81

At the end of simulation Ratio

1989 1.

3804 1.91

7256 3.65

16,555 8.32

31,447 15.81

and genotype and based on extensive experimental data in the way we have here, nor have they been based on data to the extent that our predictions are. The experimental data sets we used were collected over a long period of time, the results were confirmed in field trials and are the basis for the current rates used to control resistant insects in Australia. In addition, no previous studies of resistance in stored grain pests have used detailed individual-based models that allow this level of biological and genetic detail to be represented (Renton, 2012), although similar approaches have been employed to investigate the evolution of resistance to Bt transgenic crops in insects in the field (Storer et al., 2003), and to herbicides in weeds (Renton et al., 2011). These kind of improved models, based on better data and including more biological detail, will help us predict the evolution of phosphine resistance in R. dominica, and weigh the merits of various management options in practice for delaying or avoiding evolution of resistance as accurately as possible. As mentioned previously, our results show clearly that extending the period of fumigation is an effective strategy for decreasing the overall survival of R. dominica. These results also have important implications for management of resistance evolution, as they indicate how best to effectively control resistant genotypes, as well as susceptible genotypes. Among the nine genotypes in the two-locus model, there were five with relatively high levels of resistance: sr, hr, rs, rh and rr. All but rr beetles were killed in our simulations, meaning that the frequency of resistance alleles is always 100%. However, in reality there is likely to be immigration of non-resistant insects from outside the storage after fumigation that will dilute the frequency of resistance alleles. In fact, previous simulations have shown that the development of resistance may be suppressed by the immigration of susceptible insects from unsprayed reservoirs of infestation and thus the movement of insects within the farm or storage facility may be very important, particularly when resistance is recessive (Sinclair and Alder, 1985; Tabashnik and Croft, 1982). We argue that this suppression via immigration is also likely to more effective when resistance is caused by a combination of two genes where both genes are required to achieve significant levels of resistance. If such immigration occurs, then the fewer rr beetles surviving at the end of fumigation, the less development of resistance evolution in R. dominica would be expected to occur. Thus, the management strategy of extending exposure period of phosphine fumigation is a very important one to control or delay the resistance evolution in this pest. For example, a PH3 treatment of 0.53 mg/l (350 ppm)  7 days is often used by industry in Australia. It may be better to instead use an equi-toxic treatment with a lower concentration for a longer time, for example, “0.17 mg/l (z112 ppm)  14 days”. We will conclude with some comments about the limitations of our model and recommendations for further work. This study was based on assumptions that the fumigation concentrations within a storage facility or silo at any time are constant at their target value and thus does not account for processes of diffusion over time through the silo, sorptionedesorption into and out of

31

the grain, and leakage out of the silo. The PH3 concentration and time duration applied would in practice be enough to kill all insects if these assumptions were true and the evolution of resistance was not already well advanced, and thus the evolution of resistance in the pest would never occur. The representation used in this study is thus a somewhat idealised spatially and temporally homogenous storage facility that can be seen to represent what managers should ideally be aiming to achieve. This simplification was appropriate to allow us to address the questions we wanted to focus on in this study. However, in reality this homogeneity of dose is likely to be difficult to achieve perfectly, and both temporal and spatial variability will be inevitable to some degree. Further studies representing spatial and temporal heterogeneity of dose are needed to investigate what effects such variability will have on the short-term and long-term efficacy of different control options and the contribution of resistance alleles to further generations. The effects of emigration and immigration of insects into the storage will also be important to consider in future studies. Acknowledgement The authors would like to acknowledge the support of the Australian Government’s Cooperative Research Centres Program. We also thank Rob Emery and Yonglin Ren for their great help in provision of raw data and information about beetle life cycles and silo fumigation, and the valuable comments and suggestions of an anonymous reviewer. References Banks, H.J., 1989. Behaviour of gases in grain storages. In: Champ, B.R., Highley, E., Banks, H.J. (Eds.), Fumigation and Controlled Atmosphere Storage of Grain, Proceedings of an International Conference Held at Singapore, 14e18 February 1989, Singapore, pp. 96e107. Bond, E.J., 1984. Manual of Fumigation for Insect Control. In: FAO Plant Production and Protection Paper 54. Rome, FAO. Cassells, J.A., Darby, J.A., Green, J.R., Reuss, R., 2003. Isotherms for Australian wheat and barley varieties. In: Black, C.K., Panozzo, J.F. (Eds.), Cereals 2003, Proceedings of the 53rd Australian Cereal Chemistry Conference, 7th to 10th Sept 2003. Glenelg, South Australia, pp. 134e137. Chaudhry, M.Q., 2000. Phosphine resistance. Pesticide Outlook 11, 88e91. Collins, P.J., Daglish, G.J., Nayak, M.K., Ebert, P.R., Schlipalius, D.I., Chen, W., Pavic, H., Lambkin, T.M., Kopittke, R.A., Bridgeman, B.W., 2000. Combating resistance to phosphine in Australia. In: Donahaye, E.J., Navarro, S., Leesch, J.G. (Eds.), Proceedings of the International Conference for Controlled Atmosphere and Fumigation in Stored Products, Fresno California, USA, pp. 593e607. Collins, P.J., Daglish, G.J., Bengston, M., Lambkin, T.M., Pavic, H., 2002. Genetics of resistance to phosphine in Rhyzopertha dominica (Coleoptera: Bostrichidae). Journal of Economic Entomology 95, 862e869. Collins, P.J., Daglish, G.J., Pavic, H., Kopittke, R.A., 2005. Response of mixed-age cultures of phosphine-resistant and susceptible strains of lesser grain borer, Rhyzopertha dominica, to phosphine at a range of concentrations and exposure periods. Journal of Stored Products Research 41, 373e385. Collins, P.J., 2006. Resistance to chemical treatments in insect pests of stored grain and its management. In: Lorini, I., Bacaltchuk, B., Beckel, H., Deckers, D., Sundfeld, E., Santos, J.P., Biagi, J.D., Celaro, J.C., Faroni, L.R., Bortolini, L.O.F., Sartori, M.R., Elias, M.C., Guedes, R.N.C., Fonseca, R.G., Scussel, V.M. (Eds.), Proceedings of the 9th International Working Conference on Stored Product Protection, 15e18 October 2006, Campinas, Brazil, pp. 277e282. Daglish, G.J., 2004. Effect of exposure period on degree of dominance of phosphine resistance in adults of Rhyzopertha dominica (Coleoptera: Bostrychidae) and Sitophilus oryzae (Coleoptera: Curculionidae). Pest Management Science 60, 822e826. Emekci, M., 2010. Quo vadis the fumigants? In: Carvalho, O.M., Fields, P.G., Adler, C.S., Arthur, F.H., Athanassiou, C.G., Campbell, J.F., Fleurat-Lessard, F., Flinn, P.W., Hodges, R.J., Isikber, A.A., Navarro, S., Noyes, R.T., Riudavets, J., Sinha, K.K., Thorpe, G.R., Timlick, B.H., Trematerra, P., White, N.D.G. (Eds.), Proceedings of the 10th International Working Conference on Stored Product Protection, 27 Junee2 July 2010, Estoril, Portugal, pp. 303e313. FAO, 1975. Recommended methods for the detection and measurement of resistance of agricultural pests to pesticides: tentative method for adults of some major species of stored cereals with methyl bromide and phosphine. FAO Method No 16. FAO Plant Protection Bulletin 23, 12e25.

32

M. Shi et al. / Journal of Stored Products Research 51 (2012) 23e32

Flinn, P.W., Hagstrum, D.W., Muir, W.E., Sudayappa, K., 1992. Spatial model for simulating changes in temperature and insect population dynamics in stored grain. Environmental Entomology 21, 1351e1356. Flinn, P., Phillips, T., Hagstrum, D., Arthur, F., Throne, J., 2001. Modeling the effects of insect stage and grain temperature on phosphine-induced mortality for Rhyzopertha dominica. In: Donahaye, E.J., Navarro, S., Leesch, J.G. (Eds.), Proceedings of the International Conference for Controlled Atmosphere and Fumigation in Stored Products, Fresno California, USA, pp. 531e539. Hagstrum, D.W., Flinn, P.W., 1990. Simulations comparing insect species differences in response to wheat storage conditions and management practices. Journal of Economic Entomology 83, 2469e2475. Herron, G.A., 1990. Resistance to grain protectants and phosphine in coleopterous pests of grain stored on farms in New South Wales. Australian Journal of Entomology 29, 183e189. Longstaff, B.C., 1988. A modeling study of the effects of temperature manipulation upon the control of Sitophilus oryzae (Coleoptera: Curculionidae) by insecticide. Journal of Applied Ecology 25, 163e175. Mills, K.A., Athie, I., 1999. The development of a same-day test for the detection of resistance to phosphine in Sitophilus oryzae (L.) and Oryzaephilus surinamensis (L.) and findings on the genetics of the resistance related to a strategy to prevent its increase. In: Zuxun, Jin, Yongsheng, Liang, Xianchang, Tan, Lianghua, Guan (Eds.), Proceedings of the 7th International Working Conference on Stored-Product Protection, Beijing, China, 14e19 October 1998, pp. 594e602. National Research Council, Committee on Strategies for the Management of Pesticide Resistant Pest Populations, 1986. Pesticide Resistance: Strategies and Tactics for Management. National Academy Press, Washington, D.C. Nayak, M.K., Collins, P.J., Pavic, H., Kopittke, R.A., 2003. Inhibition of egg development by phosphine in the cosmopolitan pest of stored products Liposcelis bostrychophila (Psocoptera: Liposcelididae). Pest Management Science 59, 1191e1196. Nayak, M.K., Holloway, J., Pavic, H., Head, M., Reid, R., Collins, P.J., 2010. Developing strategies to manage highly phosphine resistant populations of flat grain beetles in large bulk storages in Australia. In: Cavalho, M.O., Fields, P.G., Adler, C.S., Arthur, F.H., Athanassiou, C.G., Campbell, J.F., Fleurat-Lessard, F., Flinn, P.W., Hodges, R.J., Isikber, A.A., Navarro, S., Noyes, R.T., Riudavets, J., Sinha, K.K., Thorpe, G.R., Timlik, B.H., Trematerra, P., White, N.D.G. (Eds.), Proceedings of the 10th International Working Conference on Stored Product Protection, 27 Junee2 July 2010, Estoril, Portugal, pp. 396e401. Pike, V., 1994. Laboratory assessment of the efficacy of phosphine and methyl bromide fumigation against all life stages of Liposcelis entomophilus (Enderlein). Crop Protection 13, 141e145. Pimentel, M.A.G., Faroni, L.A., Totola, M.R., Guedes, R.N.C., 2007. Phosphine resistance, respiration rate and fitness consequences in stored-product insects. Pest Management Science 63, 876e881.

Price, N.R., Bell, C.H., 1981. Structure and development of embryos of Ephestia cautella (Walker) during anoxia and phosphine treatment. International Journal of Invertebrate Reproduction 3, 17e25. Rajendran, S., 2000. Inhibition of hatching of Tribolium castaneum by phosphine. Journal of Stored Products Research 36, 101e106. Renton, M., Diggle, A., Manalil, S., Powles, S., 2011. Does cutting herbicide rates threaten the sustainability of weed management in cropping systems? Journal of Theoretical Biology 283, 14e27. Renton, M., 2012. Shifting focus from the population to the individual as a way forward in understanding, predicting and managing the complexities of evolution of resistance to pesticides. Pest Management Science (accepted 07.04.12., doi pending). Schlipalius, D.I., Cheng, Q., Reilly, P.E., Collins, P.J., Ebert, P.R., 2002. Genetic linkage analysis of the lesser grain borer Rhyzopertha dominica identifies two loci that confer high-level resistance to the fumigant phosphine. Genetics 161, 773e782. Schlipalius, D.I., Chen, W., Collins, P.J., Nguyen, T., Reilly, P.E., Ebert, P.R., 2008. Gene interactions constrain the course of evolution of phosphine resistance in the lesser grain borer, Rhyzopertha dominica. Heredity 100, 506e516. Shi, M., Renton, M., 2011. Numerical algorithms for estimation and calculation of parameters in modelling pest population dynamics and evolution of resistance in modelling pest population dynamics and evolution of resistance. Mathematical Biosciences 233, 77e89. Shi, M., Renton, M., Collins, P.J., 2012. Mortality estimation for individual-based simulations of phosphine resistance in lesser grain borer (Rhyzopertha dominica). In: Chan, F., Marinova, D., Anderssen, R.S. (Eds.), MODSIM2011, 19th International Congress on Modelling and Simulation. December 2011, Perth, Australia, pp. 352e358. http://www.mssanz.org.au/modsim2011/A3/shi.pdf. Shi, M., Renton, M., Ridsdill-Smith, J., Collins, P.J., Constructing a new individual-based model of phosphine resistance in lesser grain borer (Rhyzopertha dominica): do we need to include two loci rather than one? Pest Science. in press. Sinclair, E.R., Alder, J., 1985. Development of a computer simulation model of stored product insect populations on grain farms. Agricultural Systems 18, 95e113. Storer, N.P., Peck, S.L., Gould, F., Van Duyn, J.W., Kennedy, G.G., 2003. Spatial processes in the evolution of resistance in Helicoverpa zea (Lepidoptera: Noctuidae) to Bt transgenic corn and cotton in a mixed agroecosystem: a biologyrich stochastic simulation model. Journal of Economic Entomology 96, 156e172. Tabashnik, B.E., Croft, B.A., 1982. Managing pesticide resistance in crop-arthropod complexes: interactions between biological and operational factors. Environmental Entomology 11, 1137e1144. Tabashnik, B.E., 1989. Managing resistance with multiple pesticide tactics: theory, evidence, and recommendations. Journal of Economic Entomology 82, 1263e1269. Winks, R.G., 1985. The toxicity of phosphine to adults of Tribolium castaneum (Herbst): phosphine-induced narcosis. Journal of Stored Products Research 21, 25e29.