A cellular automaton to describe a pest invasion within the northern Andes

François Rebaudo1, Véronica Crespo-Pérez1,2, Olivier Dangles1,2

1IRD, UR 072, LEGS, UPR 9034, CNRS 91198 Gif-sur-Yvette Cedex and Université Paris-Sud 11, 91405 Orsay Cedex, France
2Pontificia Universidad Católica del Ecuador, Facultad de Ciencias Naturales y Biológicas, Laboratorio de Entomología, Quito, Ecuador



The SimPolilla model was developed to describe the invasion and diffusion of the potato tuber moths (PTM) (Tecia solanivora, Phthorimaea operculella and Symmetrischema tangolias, Gelechiidae, Lepidoptera), tiny moths that invaded the agricultural landscape of the North Andean region in the last decades. The larvae of these moths are serious pests of potatoes, one of the main food crops of the region. The model was developed and validated in a pilot region of central Ecuador but was built to be applicable to a much wider geographic range in North Andes.

The model description follows the ODD protocol for describing individual- and agent-based models (Polhill et al. 2008; Grimm et al. 2006; Grimm and Railsback 2005) and cellular automaton (Grimm et al. 2006, appendix A p136-147). Note that in the case of cellular automaton, some of the design concepts of the ODD protocol do not apply.

State variables and scales

The model is based on biological and ecological rules derived from field and laboratory experimental data for the three PTM species (Dangles & Carpio 2008; Dangles et al. 2008; Roux and Baumgärtner 1997). The Simiatug valley, used as a pilot region to build the model, is located in the province of Bolívar, in the central highlands of Ecuador. We focused on a study area of 20 x 20 km represented by a grid of 1,600 cells with a cell size of 0.25 km². Each cell i is characterized by its quality of habitat Qi i.e. the quantity of food resources available for the moth larvae. Qi is a boolean variable which depends on the land use (ie. crops or other uses such as woods or highlands). Each cell is also characterized by a range of temperature per cell i and month j Ti,j (in ºC), an average amount of precipitation per cell i and month j Li,j (in mm), and a mean elevation per cell i Ei (in m.a.s.l.).

Table 1. State and higher-level variables of SimPolilla

Name of variable Description Parameter Units
State variables
Elevation Elevation on the study zone per cell i Ei meters above sea level
Temperature Average temperature per cell i and month j Ti,j °C
Precipitation Average amount of precipitation per cell i and month j Li,j mm
Habitat quality Presence of potato cultures in cell i Qi Boolean
Higher-level variables
Moth abundance Immature abundance in cell i Ji Number of individuals
Moth abundance Adult abundance in cell i Mi Number of individuals
Moth abundance Gravid female abundance in cell i Gi Number of individuals

The higher-level entities are based on the number of insects of the three PTM species. Each time step represents one PTM generation based on T. solanivora life cycle duration (i.e. about 3 months at 15°C). An adjustment is made on the two other species so that each step corresponds to one PTM generation. The 500 per 500 m scale for cells was chosen for fitting the level of precision we have concerning PTM dispersion, based on available knowledge on moth dispersion (Cameron et al. 2002 [B]). Temperatures, precipitations and elevations have a 1 per 1 km resolution. Inside a square of four cells, these parameters have the same value.

Process overview and scheduling

In this section, we briefly describe the processes and scheduling of our model. Details are given in the submodels section. Each process is presented according to its sequence proceeding and in the order at which state variables are updated. Each time step is one T. solanivora generation.

Table 2. Processes of SimPolilla model

Process Submodel
State variable update State variable update
Stochastic temperature Stochastic temperature
Stochastic rainfall Stochastic rainfall
PTM mortality Crude mortality
Temperature dependent mortality
Precipitation dependent mortality
PTM dispersal PTM neighborhood dispersal
PTM reproduction Mating rate
Sex ratio
Temperature dependent fecundity

Process: state variables update

Each state variable corresponding to real data (Almanaque Electrónico de Ecuador by Alianza Jatun Sacha - CDC, digitised by DINAREN, 2003 ; Hijmans et al., 2005), is imported from individual files (one per layer), so that SimPolilla may be easily adapted to other regions.

Process: stochastic temperature

Mean temperature is transformed according to a stochastic factor (Box & Muller 1958).

Process: stochastic rainfall

Three consecutive monthly precipitations are randomly chosen during a given step.

Process: PTM mortality

PTM population is updated according to crude, temperature and precipitation mortality.

Process: PTM dispersal

PTM disperse through the territory from one cell to another by diffusion.

Process: PTM reproduction

PTM populations are updated according to biological rules (mating rate, sex ratio, fecundity). A correction is made over fitness on the two other PTM species to adjust time step based on T. solanivora life cycle.


In Simpolilla model, moths' implicit objective is to infest the considered landscape by maximizing dispersal speed. Emergent key results are level of infestation in each cell and infestation speed. Interspecific interactions are not taken into account in this model and a stochastic factor over temperature and rainfalls are included mimic actual climatic variation.


Initialization and input

The environment is based on the Simiatug agricultural region (Bolivar, central part of Ecuadorian Andes), with temperature, precipitation and elevation from available data with a 1 km² resolution. The variables were obtained from the Worldclim data set (Hijmans et al. 2005) and the BINU Project (Biodiversity Indicators for National Use, Ministerio del Ambiente Ecuador and EcoCiencia 2005). Both temperature and precipitation data corresponded to the means of the period 1961-1990 (Hijmans et al. 2005). The cellular automaton is a 4-connex square shaped grid, with closed boundaries as we are considering an existing geographical location. At the beginning of each simulation, PTM inoculums are placed in the Simiatug village and spread is observed and recorded for each species.


In this section we describe the submodels given in table 2.

Climatic driver of PTM dynamic

In order to feed the model with real climate variables, we chose to introduce a stochastic factor in the model (see also Sikder et al. 2006). We use the polar form of the Box-Muller transformation (Box and Muller 1958), to generate a Gaussian random number, based on climatic data from the region (see Dangles et al. 2008 - appendix A). Random number used is based on random procedure in VisualWoks (VisualWorks® NonCommercial, 7.5 of April 16, 2007. Copyright © 1999-2007 Cincom Systems, Inc. All Rights Reserved.).


As for temperature, we choose to introduce a stochastic factor to obtain a stochastic precipitation. Using a random number j from 1 to 12 (VisualWorks® NonCommercial, 7.5 of April 16, 2007. Copyright © 1999-2007 Cincom Systems, Inc. All Rights Reserved.), we use the monthly amounts of precipitation corresponding to three following months.

PTM life dynamics

Data on development and survival for immature stages (eggs, larva, and pupa) and on fecundity for adults were acquired from two sources. First we used published data from laboratory experiments performed in the Andean region (Notz et al 1995; Dangles et al. 2008). Second we used data obtained within the last 8 years at the Entomology Laboratory at the PUCE (Pollet, Barragan and Padilla, unpublished data).

Crude mortality

Overall force of mortality among a population is the sum of crude cause-specific forces. Here we consider dispersal related mortality (kd) and predation (ke) (see Roux, 1993; Roux et al, 1997) for P. operculella.

Sd,e(x) = e-(kd+ke)*x          (1)

Temperature dependent mortality

Temperature is the most basic controller in poikilothermic organisms (Zaslavski et al. 1988). Survival rate under laboratories conditions has been studied for the three PTM species, using a temperature dependent kinetic model (see figure 3 in Dangles et al., 2008).

Precipitation dependent mortality

Moth abundance was reduced by 80 %, when the cumulated rainfall during 3 consecutive months was higher than 600 mm (i.e. about 2.4 times more rainfall than on normal years).

Adult moth dispersal: neighborhood dispersal or dispersal per diffusion

The dispersal of the moths takes place during the adult stage when adults fly in order to find mates and/or suitable oviposition sites in potato fields (Barragán 2005). To include neighborhood dispersal into our model we considered the two factors:

  • 1) the density dependent nature of emigration rate (see Eizaguirre et al. 2004 for the Mediterranean corn borer [Sesamia nonagrioides, Lepidoptera: Noctuidae]; BenDor and Metcalf 2006 and BenDor et al. (2006) for the Emerald ash borer [Agrilus planipennis, Coleoptera: Buprestidae]) and

  • 2) the decrease in emigration rate with increasing distances (see Cameron et al. 2002 [a],[b] for P. operculella).

These factors were integrated into our cellular automata for the three species, through:

  • 1) Fraction of adults emigrating from cell i (VMi) as a function of adult density.- The maximum number of adults (K, the carrying capacity) in each cell of our model was set to 1000 individuals. Based on BenDor and Metcalf (2006) and BenDor et al. (2006) we assumed that the fraction of adults emigrating per generation (VMi), with respect to population density, followed an S-shaped curve, which levels out as density approaches 50 % of the carrying capacity (Fig. 1). We calculated VMi using equation 2 as follows:

    VMi = 0.5 / { 1 + exp [ - (Mi - ß) / µ ] }          (2)

    where Mi is the number of adults in cell i, ß is the rate of increase in migration with density (transition center) and µ is the transition width. Due to the absence of data for potato moth about the parameters of the S-shaped curve, we fixed arbitrarily ß = 500 and µ = 75, i.e. we assumed a symmetric pattern of the curve from a moth density from 0 to half of K. A previous sensitivity analysis revealed that these two parameters had little influence on the overall dispersion of moths (Rebaudo and Dangles 2008).

  • 2) Emigration rate (Pdist) as a function of distance.- Following Cameron et al. (2002 [a],[b]) for P. operculella we calculated the probability of moths flying a given distance (d) with equation (3):
  • Pdist = exp-e*d          (3)

    where e is a fixed parameter of emigration rate (see below, Fig. 2). As stated by Cameron et al. (2002 [b]), little is known about the flight capacity of potato tuber moths, and the little information available is contradictory, with some reports suggesting limited flight capacity (Fenemore 1988, Palacios et al. 1997), and others indicating that adults are good fliers (Yathom 1968). These discrepancies can probably be understood by the fact that, as for many insects, flight capacities, especially at long distances, are closely related to wind conditions (Reed 1971, Goldson and Emberson 1977, Foot 1979). Cameron et al. (2002 [b]) measured a maximum dispersal distance of 250 m for P. operculella, a value that we used for T. solanivora and Symmetrischema tangolias, as we assumed that wind transport should be similar between the three species. Following Cameron et al. (2002 [b]), we fixed e to 0.015. Concerning dispersal at short distance, comparative observations of flight capabilities between P. operculella, Symmetrischema tangolias and T. solanivora in Ecuador revealed that the latter is a much worst flier than the former (Mesías and Dangles, pers. obs.) Although our model may overestimate moth neighborhood dispersal we preferred this to below estimation.
Given the discrete nature of cellular automata, in both time and space, we could not include dispersal as a continuous variable into SimPolilla. Therefore, for each time step of the simulation, we calculated the probability of adult moths crossing the border of their current cell (Pcross). First, we assumed that moths moved inside the cell either horizontally or vertically and that they flew to the closest of the four neighboring cells. Then we considered that the probability of them crossing the cell depended on the distance flown and on cell's dimensions as follows:

Pcross(A,C,d) = { (A - (C-2*d)²) / A } with d[0:250]          (4)

where A represents the cell's surface and C its length. Pcross was then multiplied by the probability of moths emigrating (eq. 3) in order to obtain the actual probability of moths leaving a cell, which we named Pleaving:

Pleaving = Pcross * Pdist          (5)

Finally, the number of moths dispersing to neighboring cells (Ndisp) was calculated as follows:

Ndisp = VMi * Pleaving          (6)

Fig. 3 shows effective dispersal rate in relation to moth density and flight distance.

Figure 1. Moth adults emigrating fraction as a function of adult density (eq.2). Carrying capacity (K) was fixed to 1000 adults per cell i.

Figure 1. Moth adults emigrating fraction as a function of adult density (eq.2). Carrying capacity (K) was fixed to 1000 adults per cell i.

Figure 2. Moth adults emigrating fraction as a function of distance as given by equation 3. Maximum flight distance was set to 250 m.

Figure 2. Moth adults emigrating fraction as a function of distance as given by equation 3. Maximum flight distance was set to 250 m.

Figure 3. Effective dispersal rate considering moth density and flight distance.

Figure 3. Effective dispersal rate considering moth density and flight distance.

Mating rate

Mating rate is correlated with age, sex ratio and weigh of individuals, but also with distance from one individual to another (Makee et al. 2001; Cameron 2004). No specific studies have been made on mating rate under natural conditions, and laboratory measurements may frequently represent an overestimation of the natural situation because laboratory females have little opportunity to avoid mating (Reinhardt 2007). However, as our cells are 500 m long, and thanks to field observation, we know that pheromones works at least on a 200 m radius, and we assume that within a cell, mating rate is equal to one no matter the density.

Sex ratio

Among PTM population, sex ratio has been studied and is 1:1 (Saour 1999). After dispersal, the remaining adult population, combined with the mating rate and the sex ratio give the gravid females population.

PTM fecundity

Although energy-partitioning models have been developed to explain the shape of insect fecundity as a function of aging (Kindlmann et al. 2001), we are not aware of any mechanistic models that describes insect fecundity as a function of temperature. Several probabilistic non-linear models to fit insect fecundity across temperature have been proposed in the literature (Roux 1993; Lactin et al. 1995; Kim and Lee 2003; Bonato et al. 2007), but none of them gave us significant results with our data (r < 0.35, Fstat < 2.01). We therefore used weighted least square (WLS) regression to find the best model that fits our fecundity data across temperature. WLS regression is particularly efficient to handle regression situations in which the data points are of varying quality, i.e. the standard deviation of the random errors in the data may be not constant across all levels of the explanatory variables, which could be the case. For the three tuber moth species, the best fit was obtained with the Weibull distribution function, which has long been recognized as useful function to model insect development (Wagner et al. 1984).

The effect of temperature upon fecundity was well described by the Weibull distribution functions (r² = 0.75, 0.83, and 0.91 for T. solanivora, S. tangolias and P. operculella, respectively). Results showed marked differences among PTM species, both in terms of total numbers of eggs laid per females and optimal fecundity temperature : the highest fecundity was found in T. solanivora, with about 300 eggs/female at 19°C, followed by S. tangolias (220 eggs/female at 15°C) and P. operculella (140 eggs at 23°C).

  • BENDOR T.K, Metcalf S.S., Fontenot L.E., Sangunett B. and Hannon B., 2006. Modeling the spread of the Emerald Ash Borer. Ecological modelling 4328
  • BENDOR T.K., Metcalf S.S., 2006. The spatial dynamics of invasive species spread. System Dynamics Review Vol. 22, No. 1, 27-50.
  • BONATO O, Lurette A, Vidal C, Fargues J 2007. Modelling temperature-dependant bionomics of Bemisia tabaci (Q-biotype). Physiol Entomol 32:50-55
  • BOX, G.E.P, M.E. Muller 1958. A note on the generation of random normal deviates, Annals Math. Stat, V. 29, pp. 610-611.)
  • CAMERON P.J., Walker, G. P.; Penny, G. M.; Wigley, P. J 2002 [a]. Movement of Potato Tuberworm (Lepidoptera: Gelechiidae) within and Between Crops, and Some Comparisons with Diamondback Moth (Lepidoptera: Plutellidae). Environ. Entomol. 31(1): 65-75.[A]
  • CAMERON P.J., Walker, G. P.; Penny, G. M.; Wigley, P. J 2002 [b]. Movement of potato moth estimated by mark recapture experiments. 55th Conference Proceedings (2002) of TheNew Zealand Plant Protection Society Incorporated. [B]
  • CAMERON P.J., A.R. Wallace, V.V. Madhusudhan, P.J. Wigley, M.S. Qureshi & G.P. Walker, 2004. Mating frequency in dispersing potato tuber moth, Phthorimaea operculella, and its influence on the design of refugia to manage resistance in Bt transgenic crops. Entomologia Experimentalis et Applicata, Volume 115, Issue 2, Page 323-332.
  • DANGLES, O. & Carpio, F C, 2008. Cuando los científicos y las comunidades andinas unen sus esfuerzos para luchar en contra de una plaga invasora. Nuestra Ciencia 10: 23-25
  • DANGLES, O., Carpio, C, Barragan, A, Zeddam, J-L. and Silvain, J.-F. 2008. Temperature as a key driver of ecological sorting among invasive pest species in the tropical Andes. Ecological Applications, in press.
  • EIZAGUIRRE M, López C. and Albajes R., 2004. Dispersal capacity in the Mediterranean corn borer, Sesamia nonagrioides. Entomologia Experimentalis et Applicata 113: 25-34.
  • FENEMORE, P. G. 1988. Host-plant location and selection by adult potato moth, Phthorimaea operculella (Lepidoptera: Gelechiidae): a review. Journal of Insect Physiology 34: 175-177.
  • FOOT, M. A. 1979. Bionomics of the potato tuber moth, Phthorimaea operculella (Lepidoptera: Gelechiidae), at Pukekohe. New Zealand Journal of Zoology 6: 623-636.
  • GOLDSON, S. L., and R. M. Emberson. 1977. Suction trap studies of the potato moth Phthorimaea operculella (Zeller) and some observations on its biology. New Zealand Journal of Agricultural Research 20:519-525.
  • GRIMM V and Railsback S F, 2005. Individual-based Modeling and Ecology. Princeton University Press, Princeton, NJ.
  • GRIMM V, Berger U, Bastiansen F, Eliassen S, Ginot V, Giske J, Goss-Custard J, Grand T, Heinz S K, Huse G, Huth A, Jepsen J U, Jørgensen C, Mooij W M, Müller B, Pe'er G, Piou C, Railsback S F, Robbins A M, Robbins M M, Rossmanith E, Rüger N, Strand E, Souissi S, Stillman R A, Vabø R, Visser U and DeAngelis D L (2006). A standard protocol for describing individual-based and agent-based models. Ecological Modelling 198 (1-2), 115-126.
  • HIJMANS R.J. Spooner D.M., 2001. Geographic distribution of wild potato species. American Journal of Botany 88: 2101-2112.
  • KIM, D-S. and J-H. Lee. 2003. Oviposition model of overwintered adult Tetranychus urticae (Acarina: Tetranychidae) and mite phenology on the ground cover in apple orchards. Exp. Appl. Acarol. 31: 191-208.
  • KINDLMANN, P., A.F.G. Dixon, and I. Dostálková. 2001. Role of ageing 1091 and temperature in shaping reaction norms and fecundity functions in insects. Journal of Evolutionary Biology 14:835-840.
  • KOBORI Y, Amano H, 2003. Effect of rainfall on a population of the diamondback moth, Plutella xylostella (Lepidoptera: Plutellidae). Appl. Entomol. Zool. 38 (2): 249-253.
  • LACTIN, D. J., N. J. Holliday, D. L. Johnson, and R. Craigen. 1995. Improved rate model of temperature-dependent development by arthropods. Environ. Entomol 868-872.13.
  • MAKEE H., et Saour G., 2001. Factors Influencing Mating Success, Mating Frequency, and Fecundity in Phthorimaea operculella (Lepidoptera: Gelechiidae). Environ. Entomol. 30(1): 31-36.
  • NOTZ, A. 1995. Influencia de la temperatura sobre la 1138 biología de Tecia solanivora (Povolny) criadas en tubérculos de papa Solanum tuberosum L. Boletín Entomológico Venezolano 11:49-54.
  • PALACIOS, M., G. Sotelo, and E. Saenz. 1997. La Polilla de la Papa Tecia solanivora (Povolny). Pages 1-22 in Primer seminario taller internacional sobre manejo integrado de Tecia solanivora. International Potato Center (CIP), INIAP, Ibarra, Ecuador.
  • POLHILL J.G, Parker D, Brown D and Grimm V 2008. Using the ODD Protocol for Describing Three Agent-Based Social Simulation Models of Land-Use Change. Journal of Artificial Societies and Social Simulation vol. 11, no. 2 3
  • REBAUDO, F., and O. Dangles 2008. Developing multi-agent based models to strengthen food security in the Northern Andes: Application to invasive agricultural pests in Central Ecuador. Technical Assistance Report for the McKnight Foundation; Quito, Ecuador.
  • REED, E. M. 1971. Factors affecting the status of a virus as a control agent for the potato moth (Phthorimaea operculella (Zell.) (Lep.: Gelechiidae)). Bulletin of Entomological Research 61: 207-222.
  • REINHARDT K., Köhler G., Webb S., and Child D., 2007. Field mating rate of female meadow grasshoppers, Chorthippus parallelus, estimated from sperm counts. Ecological Entomology, Volume 32, Issue 6, Page 637-642.
  • ROUX O, 1993. Population ecology of potatoes tuber moth Phthorimaea operculella (Zeller) (Lepidoptera: Gelechiidae) and design of an integrated pest management program in Tunisia. Thesis ETH No 10120, Zurich.
  • ROUX O, Baumgärtner J, 1997. Evaluation of mortality factors and risk analysis for the design of an integrated pest management system. Ecological modeling 109, 61-75.
  • SAOUR G., 1999. Susceptibility of potato plants grown from tubers irradiated with stimulation doses of gamma irradiation to potato tuber moth, Phthorimaea operculella Zeller (Lep., Gelechiidae). Journal of Applied Entomology, Volume 123, Issue 3, Page 159-164.
  • SIKDER I.U. Sikder, Mal-Sarkar S, Mal T.K., 2006. Knowledge-Based Risk Assessment Under Uncertainty for Species Invasion. Risk Analysis, Vol. 26, No. 1.
  • TANHUANPÄÄ M, Ruohomaki, K., Kaitaniemi, P. 2003. Influence of adult and egg predation on the reproductive success of Epirrita autumnata (Lepidoptera: Geometridae). Oikos 102: 263-272.
  • WAGNER, T.L.; WU, H-I; Sharpe, P.J.H.; Schoolfield, R.M.; Coulson, R.N., 1984. Modeling insect development rates: a literature review and application of a biophysical model. Annals of the Entomological Society of America, Volume 77, Number 2 , pp. 208-225(18).
  • WAKISAKA S., Tsukuda R and Nakasuji F, 1989. Effects of Natural Enemies, Rainfall, Temperature, Host Plants on Survival and Reproduction of Diamondback Moth. Otsuka Chemical Co., Ltd., Naruto Research Center, Naruto Tokushima 772, Japan and Laboratory of Applied Entomology, Faculty of Agriculture, Okayama University, Tsushimanaka, Okayama 700, Japan.
  • YATHOM, S. 1968. Phenology of the tuber moth, Gnorimoschema operculella Zell., in Israel in the spring. Israel Journal of Agricultural Research 18:89-90.
  • ZASLAVSKI V. A. 1988. Insect development: photoperiodic and temperature control. Springer, Berlin.

Modeling Notebook: code testing and sensitivity analysis for SimPolilla© model

APPENDIX A: code testing

APPENDIX B: sensitivity analysis

Corresponding author : François Rebaudo


Le Cirad Centre de coopération internationale en recherche agronomique pour le développement
Informations légales © Copyright Cirad 2001-2015