M ODELING THE SPREAD AND GEOGRAPHIC DISTRIBUTION OF INVASIVE TERMITES IN FLORIDA By FRANCESCO TONINI A DISSERTATION PRESENTED TO THE GRADUATE SCHOOL OF THE UNIVERSITY OF FLORIDA IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY UNIVERSITY OF FLORIDA 2014
2014 Francesco Tonini
To my family and friends
4 ACKNOWLEDGMENTS I would like to thank my colleagues and friends for supporting me throughout this long journey. A special thank goes to Dennis Zielstra, Adam Benjamin, Majid Alivand, and Sreten Cvetojevic I would also like to thank my family and close friends for their constant support and for believing in me. I am thankful to my advisor, Dr. Hochmair, for giving me the opportunity of pursuing my PhD and for his advice during these years.
5 TABLE OF CONTENT S page ACKNOWLEDGMENTS ................................ ................................ ................................ .. 4 LIST OF TABLES ................................ ................................ ................................ ............ 7 LIST OF FIGURES ................................ ................................ ................................ .......... 8 LIST OF ABBREVIATIONS ................................ ................................ ........................... 10 ABSTRACT ................................ ................................ ................................ ................... 11 CHAPTER 1 INTRODUCTION ................................ ................................ ................................ .... 13 Objectives ................................ ................................ ................................ ............... 15 Significance of the Study ................................ ................................ ........................ 16 Dissertation Outline ................................ ................................ ................................ 16 2 SIMULATING THE SPREA D OF AN INVASIVE TER MITE IN AN URBAN ENVIRONMENT USING A STOCHASTIC INDIVIDUA L BASED MODEL .............. 21 Study Background ................................ ................................ ................................ .. 21 Materials and Methods ................................ ................................ ............................ 25 Model Design ................................ ................................ ................................ ... 25 Purpose of model ................................ ................................ ....................... 25 Entities, state variables and scales ................................ ............................ 25 Process overview and scheduling ................................ .............................. 26 Design concepts ................................ ................................ ........................ 27 Initialization ................................ ................................ ................................ 27 Input data ................................ ................................ ................................ ... 28 Submodels ................................ ................................ ................................ 28 Sensitivity Analysis ................................ ................................ ........................... 31 Model Parameterization ................................ ................................ ................... 32 Application of Model to Specific Study Area and Data ................................ ............ 36 Resu lts and Discussion ................................ ................................ ........................... 37 3 STOCHASTIC SPREAD MO DELS: A COMPARISON B ETWEEN AN INDIVIDUAL BASED AND A CELL BASED MODEL FOR ASSE SSING THE EXPANSION OF INVASIV E TERMITES OVER A LA NDSCAPE ........................... 51 The Integration of GIS with Spread Models ................................ ............................ 51 Materials and Methods ................................ ................................ ............................ 53 Model Design ................................ ................................ ................................ ... 54 Framework and Modules ................................ ................................ .................. 59
6 Results ................................ ................................ ................................ .................... 60 Homogeneous Landscape ................................ ................................ ................ 60 Case Study: Dania Beach, Florida ................................ ................................ ... 62 Comparison of Model Results and Discussion ................................ ........................ 62 Advantages and Limitations of the Model ................................ ............................... 65 4 PREDICTING THE GEOGR APHICAL DISTRIBUTION OF TWO INVASIVE TERMITE SPECIES FRO M OCCURRENCE DATA ................................ ............... 77 Species Distribution Models ................................ ................................ .................... 77 Materials and Methods ................................ ................................ ............................ 79 Study Area and Species Data ................................ ................................ .......... 79 Predictor Variables ................................ ................................ ........................... 80 Future Scenarios ................................ ................................ .............................. 81 Modeling Approaches ................................ ................................ ....................... 83 Maximum entr opy approach ................................ ................................ ....... 83 Bayesian for presence only data (BPOD) approach ................................ .. 85 Evaluation of model performance ................................ .............................. 87 Sampling scheme ................................ ................................ ...................... 87 Results ................................ ................................ ................................ .................... 88 Discussion ................................ ................................ ................................ .............. 90 5 CONCLUSIONS ................................ ................................ ................................ ... 101 APPENDIX A SENSITIVITY ANALYSIS ................................ ................................ ..................... 104 B SUPPLEMENTARY DATA ................................ ................................ .................... 106 C COMPARIS ON OF MODEL OUTCOMES ................................ ............................ 111 D PRINCIPAL COMPONENTS LOADINGS ................................ ............................. 112 E AST MODEL COEFFICIENTS ................................ ................................ .............. 113 F FST MODEL COEFFICIENTS ................................ ................................ .............. 114 LIST OF REFERENCES ................................ ................................ ............................. 115 BIOGRAPHICAL SKETCH ................................ ................................ .......................... 123
7 LIST OF TABLES Table page 2 1 Model parameters: abbreviations, definitions, and their baseline values. ........... 50 3 1 Model ecological parameters used in both modeling approaches and their values. ................................ ................................ ................................ ................ 76 4 1 Environmental variables used and their data source for both AST and FST. ..... 98 4 2 List of main models used for AST, their sampling sensitivity, and information criteria. ................................ ................................ ................................ ................ 99 4 3 List of main models used for FST, their sampling sensitivity, and information criteria. ................................ ................................ ................................ .............. 100 A 1 Sensitivity analysis results extracted for some years of the simulation. ............ 104 C 1 Simulation results over a homogeneous landscape. ................................ ........ 111 D 1 Environmental and climat e variables used in the study and principal component analysis (PCA) loadings for all six extracted main axes, accounting for about 80% of the total variation in the original dataset. ............. 112 E 1 .......... 113 F 1 .......... 114
8 LIST OF FIGURES Figure page 1 1 Individuals of Coptotermes formosanus Shiraki. ................................ ................. 18 1 2 Individuals of Coptotermes gestroi (Wasmann). ................................ ................. 19 1 3 Individuals of Nasutitermes corniger (Motschulsky). ................................ ........... 20 2 1 Structure of the initialization steps involved in the termite dispersal simulation model. ................................ ................................ ................................ ................. 43 2 2 Core subprocesses involved in the individual based simulation algorithm at any generic time step. ................................ ................................ ........................ 44 2 3 Location of samples of N corniger collected during a field survey in 2003 ........ 45 2 4 Snapshot of the areas predicted as infested with N. corniger by the baseline dispersal simulation model. ................................ ................................ ................ 46 2 5 Sensitivity analysis charts for the termite dispersal simulation model. ................ 47 2 6 Model sensitivity to the SURV parameter. ................................ .......................... 48 2 7 Evaluation of the termite dispersal simulation model.. ................................ ........ 49 3 1 Schematic example of termite reproductives and the pheromone attraction radius within which a heterosexual individual must b e found to mate after a dispersal flight. ................................ ................................ ................................ .... 67 3 2 A sampling of empirical histograms for the number of new termite colonies formed b y N reproductives (alates). ................................ ................................ .... 68 3 3 Conceptual diagram of the main modules used in the simulation algorithm at any generic time step. ................................ ................................ ........................ 69 3 4 Schematic representation of the focal ("moving window") operation. ................. 70 3 5 Simulation outcome comparison between cell based and individual based models.. ................................ ................................ ................................ .............. 71 3 6 Predicted areas (in hectares) of infestation for the cell based model (solid line) and the individual based model (dashed line) for different values of N randomly distributed initial colonies. ................................ ................................ ... 72 3 7 Computational time (in hours) of the cell based model (solid line) and the individual based model (dashed line) for different values of N randomly distributed initial colonies. ................................ ................................ ................... 73
9 3 8 Snapshot of the areas predicted as infested with N. corniger by the simulation m odels with sampled termite locations in 2003 and the period 2004 2012. ................................ ................................ ................................ ......... 74 3 9 Habitat representation in the IBM and the cell based m odel. ............................. 75 4 1 Florida occurrences of AST (green) and FST (purple). ................................ ....... 94 4 2 Current predicted probabilities of presence for AST. ................................ .......... 95 4 3 Current predicted probabilities of presence for FST. ................................ .......... 96 4 4 Current and average projected probabilities of presence for the 2050s time period. ................................ ................................ ................................ ................ 97 B 1 Empirical histograms for the number of new colonies formed by N reproductives (alates) ................................ ................................ ...................... 106
10 LIST OF ABBREVIATIONS AIC Akaike Information Criterion AICc Akaike Information Criterion Corrected AUC A rea U nder the R eceiver O perating C haracteristic C urve AST Asian Subterranean Termite BPOD Bayesian for P resence O nly D ata CA Cellular Automata CCAFS Climate Change Agriculture and Food Security FST Formosan Subterranean Termite GCM Atmospheric Oceanic Global Circulation Models IBM Individual Based Model IPCC Intergovernmental Panel on Climate Change Maxent Maximum Entropy NLCD National Land Cover Database PCA Principal Components Analysis SDM Species Distribution Model SEPM S patially E xplicit P opulation M odels
11 Abstract of Dissertation Presented to the Graduate School of the University of Florida in Partial Fulfillment of the Requirements for the Degree of Doctor of Philosophy M ODELING THE SPREAD AND GEOGRAPHIC DISTRIBUTION OF INVASIVE TERMITES IN FLORIDA By Francesco Tonini May 2014 Chair: Hartwig Henry Hochmair Major: Forest Resources and Conservation Invasive termites are destructive insect pests that cause billions of dollars in property damage every year. To be able to mitigate the damaging effects of termite infestation this research investigates the spread and geographic distribution of invasive termites in Florida. In order to assess the ir spatiotemporal spread over a landscape t wo s patially explicit stochastic simulation models were developed : an individual based model (IBM) and a cell based model. The setting of a case study of an invasive termite, Nasutitermes corniger (Motschulsky), was used to simulate the spread of the species in Dania Beach, FL The Monte Carlo simulation technique was us ed in both models to assess outcome uncertainty A sensitivity analysis was carried out for the IBM using a set of model realizations describing potential areas of infestation in order to assess the importance of each ecological parameter included in the m odel Results showed that the areas with the highest predicted likelihood of infestation in both simulation models matched well the empirical locations discovered over the years within the chosen time window.
12 Frequently, the only data available on a speci es is the location of its occurrence (presence only data) In order to predict the potential habitat of invasive termite species under both current and future climate change scenarios several statistical models were considered usin g two approaches for presence only data : a recently developed Bayesian linear logistic regression approach adjusted for presence only data and the widely used maximum entropy approach (Maxent) A case study of two destructive invasive termite species in Florida Coptotermes ge stroi (Wasmann) and Coptotermes formosanus Shiraki, is presented to predict the ir current and future biotic distributions by using their occurrence records Results showed that the predicted distributions of both C gestroi and C formosanus are strongly l inked to urban development and that f uture climate warming seems to affect their projected probability of presence to a lesser extent than population growth All models presented in this study can help state and local regulatory agencies to understand the response of invasive termites to different environmental conditions and anticipate the ir spread by develop ing successful early detection, quarantine, or eradication plans based on the predicted areas of infestation
13 CHAPTER 1 INTRODUCTION Termites are serious insect pests of the urban environment and are responsible for caus ing billions of dollars worth of damage to structural lumber ( Edwards and Mill, 1986 ) Recently, an estimate of the worldwide annual control and repair cost was found to be around $40 billion ( Rust and Su, 2012 ) Human transport is responsible for spreading destructive exotic termites from their native range to far reaching pa rts of the globe. Worldwide, over a dozen termite species have become established across oceans via colonies onboard infested boats or by the movement of infested goods. In Florida no less than six non endemic termite species have become established in par ts of the state ( Scheffrahn et al., 2002 ) Once established, the natural termite dispersal of an exotic species proceeds slowly. Termite colonies require many years to mature, and once the first crop of winged reproductives leaves the colony, they are unable to fly more than a few hundred meters fr om the parent colony ( Evans, 2011 ) Anthropogenic dispersal, i.e. human assisted, is far more rapid and progresses at cruising or highway speeds. The latter is far more complex a nd unpredictable than the natural dispersal, since occurrences have been recorded at localities of anthropogenic introduction beyond the endemic range of a given termite species. However, only if the climatic conditions are suitable, will the introduced sp ecies flourish. Termites are not currently found across most of the extreme northern part of the United States and Canada However, some published articles state that it will probably not take long time until termites will start establishing in those areas ( Peterson, 2010 ) Such a quick establishment would have to rely on anthropogenic movement and rapid climate change
14 The University of Florida Termite Collection at the Fort Lauderdale Research and Education Center is the largest and newest repository of exotic an d native termites in the world. Over 32,000 geo referenced and identified colony samples are included in the database The majority of the recorded samples are located in North America and the Caribbean. In this work, a special focus w as given to three in vasive destructive termite species in Florida: the Formosan subterranean termite (FST), Coptotermes formosanu s Shiraki (see Figure 1 1) the Asian subterranean termite (AST), Coptotermes gestroi (Wasmann) (see Figure 1 2) and Nasutitermes corniger (Motsch ulsky) (see Figure 1 3) It has been suspected that FST and AST have been introduced to and dispersed throughout Florida by sailboats and yachts ( Hochmair and Scheffrahn, 2010 ) FST has been introduced in the continental United States during the 1960s and i n Florida it is found in most urban centers statewide ( Su and Scheffrahn, 2000 ; Su and Scheffrahn, 1997 ) AST has only recently been introduced in the continental United States and has established only in south east Florida ( Scheffrahn and Su, 2005 ) N cornige r is a nonendemic species in the United States and its discovery in Dania Beach, Fl orida, represents the first record of a land based establishment by a higher termite (Family Termitidae) in the continent ( Scheffrahn et al., 2002 ) The S tate of Florida was chosen as a common study region for all the aforementioned termite species given the following reasons: (i) relatively homogeneous landscape, hence ideal conditions for developing and testing spatially explicit simulation models; (ii) large climate variability along the north south direction, ranging from temperate continental in the panhandle to oceanic subtropical in the Keys ; (iii) rare
15 simultaneous pres ence of AST and FST which, at present, occur sympatrically only in south Florida, Hawaii, and Taiwan; and (iv) large amount of geo referenced records available from the University of Florida Termite Collection Objectives The main goal of this study was to investigate the spread and geographic distribution of invasive termites in order to mitigate the damaging effects of an infestation and better understand their response to both environmental and climatic conditions A central h ypothesis is that the natura l range expansion of termite pests will proceed extremely slowly over time, taking hundred/thousand years to cover a few hundred miles. Moreover, the spread will probably occur at a much slower speed than it takes for suitable habitat boundaries to change under climate change scenarios. The anthropogenic dispersal will not be considered for the purpose of this work since it is practically an unpredictable event. It is also hypothesized that shape and boundaries of suitable habitats affect the dispersal spee d. For example, street networks, canals accessible via boat, lakes, patches of open grassland, etc. might have an influence on the spread rate and flight direction The aforementioned hypotheses are tested herein by pursuing three main objectives: 1. Study th e potential speed and dispersal rate at which a n introduced termite pest might travel by natural means over a landscape 2. Predict the potential habitat of C. formosanus and C gestroi over Florida under both current and future climatic conditions. 3. Develop models that can be applied to other termite species within the constraints of their own specific natural histories
16 Significance of the Study This study tries to address important questions on the spread and distribution of invasive termite species, such as : w here could they ultimately become sustained pests and how fast, once established, can they continue to disperse on their own? How would they change and adjust their natural geographic r ange under different climate change scenarios, i.e. where would they re establish, and how fast would they spread to new geographic regions? Models predicting the areas of infestation following initial introduction of an invasive species are important for an authoritative response and early detection is critical for successful quarantine or eradication efforts. As for most organisms, the distribution of exotic termites is limited by climate and habitat. No exotic termite has yet saturated its non native ran ge except for a few small island situations, e.g., Barbados 1 Predictions for geographic establishment and spread of invasive termite pests can be used by government regulators and the pest control industry to either anticipate their arrival or mitigate th e devastating effect of an infestation in terms of economical costs and damage caused to buildings Dissertation Outline The c hapters of this dissertation are self contained journal articles The first objective of the dissertation was accomplished by deve loping two stochastic simulation models. An individual based model is described and tested in C hapter 2, while a cell based model is described, tested, and compared to the latter model in C hapter 3. The second objective of the dissertation was accomplished by considering several advanced statistical models using two main approaches for presence only data Chapter 4 1 Scheffrahn, personal communication.
17 describes the aforementioned statistical models and considers both historical climatic conditions and a set of future scenarios to be tested. Th e third objective of the dissertation was accomplished by developing all models and approaches in a way that they can easily be generalized to other termite species after considering the constraints of their own specific natural histories. Finally, conclu sions on the advantages limitations and implications of th is work are provided in C hapter 5.
18 Figure 1 1 Individuals of Coptotermes formosanus Shiraki. A) W inged reproductive (alate). B) S oldier. C) W orker.
19 Figure 1 2 Individuals of Coptotermes gestroi (Wasmann). A) W inged reproductive (alate). B) S oldier. C) W orker.
20 Figure 1 3 Individuals of Nasutitermes corniger ( Motschulsky ) A) S oldier. B) W orker (major morph) C) W orker (minor morph).
21 CHAPTER 2 SIMULATING THE SPREA D OF AN INVASIVE TER MITE IN AN URBAN ENVIRONMENT USING A STOCHASTIC INDIVIDUA L BASED MODEL Study Background The primary anthropogenic means by which termites are transported between continents and islands is by maritime vessel ( Scheffrahn and Crowe, 2011 ) Over a dozen exotic termite species have become established worldwide ( Evans, 2011 ) of which six can be found in Florida ( Scheffrahn et al., 2002 ) Termites are destructive insect pests that cause billions of dollars in pr operty damage every year ( Edwards and Mill, 1986 ) Once a species is established, th e natural dispersal of termite colonies proceeds slowly. Termite colonies typically require 4 6 years to mature, and once the first group of alates (winged reproductives) leaves the colony, they are unable to fly more than a few hundred meters from the par ent colony ( Husseneder et al., 2006 ; Messenger and Mullins, 200 5 ; Mill, 1983 ) Anthropogenic or long distance movements lack predictabili ty. Specifically, the nesting core of a termite colony (reproductives, brood, and most foragers) must be moved intact and both a water and food source must be associated with the core during movement ( Hochmair and Scheffrahn, 2010 ) The inherent complexity of a physical environment limits the applicability of mathematical models for realistic dispersal modeling of invasive species, and practical predictions are difficult to obtain ( Pitt, 2008 ) Analytical methods commonly used to model dispersal in the past and in some cases up to the present include: (i) simple reaction diffusion mo dels ( Fisher, 1937 ) which ignore any spatial interaction between individuals and do not consider single dispersal events; (ii) mixed diffusion population
22 growth models, which include a per capita growth parameter ( Okubo, 1980 ; Skellam, 1951 ) or several demographic variables ( van den Bosch et al., 1990 ) ; and (iii) integro differenti al equation models, which separate population dynamics and dispersal into two stages ( Neubert et al., 1995 ) More recently, computer intensive approaches, such as spatially explici t population models (SEPMs), have been able to incorporate both ecological/biological information at a population level with underlying habitat differences ( Wiegand et al., 2004 ) Computer simulations seek to imitate the dynamics of various real world processes ( Steyaert, 1993 ) rather than solving sets of equations. Simulation models are either deterministic or stochastic. The first model type gives a fixed output for a given set of input data and model parameters while the second model type includes at least o ne stochastic process and provides a probabilistic outcome ( Law and Kelton, 1982 ) The i ntrinsic dynamic component of a computer simulation provides the ability to estimate the rate at which an invading species is likely to occupy suitable areas. However, such models may represent a poor choice in cases where established populations are restr icted to distinct areas of suitable habitat, since assuming universal dispersal abilities may not reflect the ability of a species to move from a current location to another potentially suitable habitat ( Peterson et al., 2002 ) Whereas simulating the spread of invasive species beyond a decade into the future may decrease the reliability of the model outcome ( Pitt et al., 2011 ) it should be noted that the invasive plant used by Pitt et al. (2011) has a much faster dispersal capacity compared to termites. Individual based models (IBMs) are able to incorporate several rules describing the interactions between individual units considering each one of them as different, both
23 physiologically and behaviorally ( Huston et al., 1988 ) The complexity of the rules increases with the total number of parameters involved in describing them. However, complexity often comes at the expense of generality, which makes it necessary to select the most appropriate modeling approach on a case by case bas is. Small spatial scales, such as urban environments, are particularly suitable for the development of IBMs, because they are complex enough to require simulation but not so large as to be unmanageable for an IBM. Also, IBMs are able to represent individua ls explicitly and incorporate biologically relevant rules that have a strong influence on the dynamics of an invasion ( Pitt, 2008 ) In this work a computer simulation was developed using a spatially explicit stochastic individual based modeling approach and use hindcasting in order to predict which areas would have been infested by an arboreal invasive t ermite, Nasutitermes corniger (Motschulsky), had no eradication plan been implemented at a particular location, Dania Beach, FL. The methodology presented herein is appropriate for more general application, such as predicting the future geographical spread or studying a different termite species after appropriate adjustments in the model parameters. Individual based simulations consider the individual organism to be a logical basic unit for modeling ecological phenomena ( Grimm and Railsback, 2005 ) E ach model was run from 2003, the year in which a first complete survey of infested locations had been conducted over the study area, until 2012. The model outcome is the predicted areas of infestation at any time step, indicating the spatial extent and dynamic evolution of the invasion. Beginning in 2003, local authorities have been trying to eradicate this pest from the original survey area. However, between 2006 and 2011, extended survey
24 procedures had to be stopped due to discontinued funds. A new recent survey conducted in 2012 found newly infested locations in areas not spotted originally and therefore not included in the eradication plan. S tate or local regulatory agencies may benefit from a model that predicts the rate and direction of termite dispersal, as it would as sist them in targeting specific areas for survey, eradication, or quarantine efforts. In the literature, only two computer simulation models have been applied to a termite species: one has been developed to determine per capita wood consumption rates of t ermite workers ( Morales Ramos and Rojas, 2005 ) while the other explored termite foraging behavior underneath the soil ( Lee et al., 2008 ) To date, no computer simulation models have been published that investigate the geographic spread of a termite infestation from a set of surveyed locations. Unlike some other recently developed spatial simulation models found in the literature for other insects ( Carrasco et al., 2010 ; Pitt, 2009 ) the human mediated dispersal component is not included because of its unpredictability and lack of calibration data. Although samples collected o ver the past 10 years do not reflect the true (i.e. natural) expansion of the species, and were collected mainly for the purpose of verifying the success of the eradication effort, it is nevertheless possible to use the newly infested locations (2012) to g round truth our simulation model. T he parameters and methods used to develop the computer simulation are described herein Results are presented together with a discussion on the relative importance of each biological parameter included in the model, follo wed by conclusions.
25 Materials and Methods Model D esign The simulation algorithm is implemented using a set of R functions ( R Development Core Team, 2013 ) and free source code is provided The model d escription follows the ODD (Overview, Design concep ts, Details) protocol ( Grimm et al., 2006 ; Grimm et al., 2 010 ) in order to make the model's logic as clear as possible. Purpose of m odel A spatial, stochastic computer simulation was developed with the purpose of gaining a deeper understanding of the rate and direction of a termite invasion by natural means over a realistic landscape, such as an urban environment. In this study, the model is also used to determine how a new invasive species in South Florida, N. corniger could have expanded from a set of surveyed locations up to the present, if no eradication plan had been implemented throughout the years. The developed simulation model may assist state or local regulatory agencies in targeting specific areas for survey, eradication, or quarantine efforts Entities, s tate v ariables and s cales The basic entities of the model are individual termite alates (dispersing propagules) and all the individual colonies they are generated from. Both alates and colonies are characterized by their continuous spat ial location specified in a Cartesian plane coordinate system. Alates are also characterized by their sex (M F), and colonies by their age (in years). A reference spatial grid was used to represent the distribution of all areas occupied by one or more term ite colonies at each time step. The grid is set to an extent of 10 km x 10 km over the urban area of Dania Beach, FL, with a resolution of 100 x 100 meters. T he chosen resolution is suitable for a few reasons such as the
26 uncertainty associated with the pre cise locations of surveyed colonies/individuals, the ( Collins, 1981 ) and because it is a suitable scale of surveillance and pest control management. In order for the simulation to be more realistic, the local urban environment was considered and the areas that are unsuitable for the establishment of a new colony were excluded such as roads, highways, non wooded fields, and water bodies. Each area wit h wood sources (e.g., buildings, trees, boats, debris, etc.) has potential for colonization. F or the chosen temporal resolution (10 years) the choice of a static habitat suitability layer should not introduce any relevant bias in the results. However, shou ld the model be run over a much longer time span, a different strategy may be considered The temporal scale is discrete and one time step represents 1 year. The model is run from 2003 (year of the first complete survey of infested locations) to 2012. Proc ess overview and scheduling Dispersal of alates is the key process in the spread of colonies, and the dispersal was simulated as a single annual event. The consequence may be an increased chance for alates to find a mate and form a new colony. However, thi s represents a necessary simplification, since typical termite dispersal is formed by a major exodus that may be preceded and/or followed by smaller flights, of unknown magnitude and timing. Many termite species initiate dispersal flights in the early stag es of the wet season and are triggered by environmental factors ( Jones et al., 1988 ; Martius, 2003 ; Nutting, 1969 ) Dispersal flights are the only means by which new colonies can form beyond the foraging territory of the mother colony. Although the model simplifies the temporal scale of the real phenomenon, single massive dispersal flights are common because: (i)
27 alates are less vulnerable as prey, as they can overwhelm predato rs by large numbers ( Nutting, 1969 ) ; and (ii) there are higher odds of finding and choosing a mate. Design c oncepts Sensing D ispersing alates (reproductives) can sense and respond t o pheromones in order to find potential mates of the opposite sex that have dispersed by chance to the nearby sites. Interaction Male and female alates interact to form new colonies. Stochasticity Both distance and direction of dispersal by alates are de termined stochastically from a probability distribution (see paragraph on Colony a ging and a late p roduction ). The sex (male or female) of a particular alate is random. Collectives Individual alates are followed during dispersal, but after a colony is form ed by two alates of the opposite sex, the colony is followed as whole rather than at the resolution of individuals. Initialization Fig. 2 1 shows a schematic representation of the steps involved for the model initialization. At the initial state, i.e. time t=1, the spatial locations of all surveyed termite colonies are stored in a dataset and assigned a random age between 0 and the maximum lifespan decided by the user. The initialization process is the same in all simulation runs. Surveyed colonies can be i mported from an external data file containing their geographic coordinates, e.g. recorded with a GPS device. In most cases, the collected samples do not identify different termite colonies, as they are taken opportunistically with the goal of spotting an i nfestation. Therefore, different termite locations may or may not belong to the same colony.
28 Input data Table 2 1 shows seven main parameters of the implemented dispersal model and their baseline values, i.e. the values assigned for the baseline simulation, which are based either on related literature findings (see subsection on Model Parameterization ) or the opinio n of termite experts. More specific information for the particular location modeled, Dania Beach, FL, is described in Application of Model to Speci fi c Study Area and Data section Submodels The simulation algorithm is composed of several modules ordered in a sequential manner and imitates the steps taken by a group of individual alates from the dispersal to the inception of a new colony. Figure 2 2 shows an overview of the main subprocesses and steps involved in the simulation at any generic time step. Each subprocess is discussed in detail below. Habitat suitability T he habitat suitability submodel checks the suitability of the underlying environment for all termite individuals after dispersal. If an individual alate falls within an unsuitable habitat, as defined by th e user, then it is eliminated. In order to include the local landscape in the simulation model and identify areas unsuitable for the establishment of new termite colonies, the following vector type spatial layers were combined in a GIS (using was obtained from the Broward County GIS data online repository ( Lelis, 2006 ) The 2011 NAVTEQ NAVSTREETS street data were used for the street network layer and created a 10 m buffe r around each line segment to model the approximate coverage of roads. Further, the Fort Lauderdale Airport area and its runways were extracted from the 2010 TomTom (formerly TeleAtlas) Multinet Dataset.
29 Because N. corniger like other invasive termite spe cies, needs a wood source as food, all agricultural field polygons were added to the collection of unsuitable areas. These polygons were extracted from a 2010 land use layer at the parcel level, which was obtained from the University of Florida GeoPlan Cen ter. The original land us e layer was compiled by the Florida Department of Transportation and contains 99 land use classes which have been collapsed into 15 classes by the GeoPlan Center ( University of Florida GeoPlan Center, 2010 ) Using the union overlay operation in ArcMap, all the GIS layers listed above were combined into a single layer denoting unsuitable habitat areas in which colonies are not able to survive. Alate dispersal The dynamics and speed of termite dispersal by natural means are controlled by several behavioral characteristics affecting the successful creation of new colonies. S uch behavioral characteristics were identified and included in t he form of model parameters to better simulate the real process. A new colony begins with a male female (i.e. king and queen) couple of unwinged alates building and sealing the nuptial chamber in a proper substrate, usually soil or wood. After a termite co lony matures, which requires about 4 years, alates are produced. All alates change their behavior in response to: (i) changes of habitat, i.e. they may fall into an unsuitable patch of land and therefore are not able to find a location to form an incipient colony; (ii) their proximity to a heterosexual mate. Alates do not adjust their behavior over time as a consequence of their experience, since they only serve the purpose of expanding the colony with a one time flight after which they either die or find a mate and become the king/queen of a new colony. Although they have eyes, alates are probably not able to predict which location will be suitable once in flight. Dispersal flights typically occur at dusk or at night
30 after a light rain and during calm weath er conditions. It is known that alates are attracted by lights, as found in mark recapture studies ( Messenger and Mullins, 2005 ) Sex pheromones have two main roles: a close range attraction and contact attraction. The former is used to unite sexual partners, the lat ter is used to maintain the contact during the tandem behavior ( Nutting, 1969 ) Alates do not release pheromones during the flight and therefore their flight behavior is not influenc ed by it. The processes that are modeled assuming they are stochastic, i.e., random, are the flight distance, flight direction, and the sex of each individual. The model output is used to spot which areas have been occupied and how often throughout all 100 runs. Colony formation The colony formation subprocess loops through each grid cell that is occupied by at least two individuals and, subsequently, through each individual. This process is necessary to check if a reproductive is able to find a heterosexu al neighbor and form a nuptial pair, where the neighborhood is defined by a circular buffer with the pre set pheromone attraction radius. If two candidate alates are matched, a new colony is created and assigned spatial coordinates of the midpoint between the two individuals. The process stops for a particular grid cell as soon as the maximum density of colonies per hectare is reached. At the end of the present subprocess, if one or more pairs of individuals are matched, new colonies are created and their s patial location is saved. Colony aging and alate production Each time step, the age of every existing colony is increased by one (aging submodel). If this value exceeds the maximum lifespan defined by the user, then the colony is eliminated from the map. After the aging subprocess, alates are generated by each existing colony based on colony age
31 (dispersal subprocess). Older colonies generate more alates, which increase the overall chances for an individual reproductive to find a mate, a suitable nesting s ite, and a location farther away from the mother colony. The dispersal subprocess also executes the following: (i) random creation of male and female individuals by sampling from a Binomial distribution, Bin(n,p), where n is the number of alates to be gene rated and p is the probability of drawing a male, (ii) random sampling of flight directions (in radians) rate), and (iv) calculation of new spatial locations X and Y (Easting and Northing) of alates derived from basic trigonometric equations, using both the simulated flight direction and flight distances. Updating the distribution of colon ies on the landscape The final subprocess (stacking colonies subprocess) stacks all new colonies created during the previous process (colony formation) with the existing ones in a dataset. Before moving to the next time step, all colonies can be saved to an external shapefile as points and further converted to a geo referenced raster grid. The raster grid allows us to overlay modeling results from multiple simulations and to compute a final occupancy envelope. At the following time step, new alates are gen erated which fly out from all mature colonies, i.e., colonies old enough to produce alates. Sensitivity A nalysis A range of values for each model parameter (Appendix A, Table A 1) was chosen to run a sensitivity analysis and assess the contribution of each one of them to the model outcome. Preliminary model runs were used to identify the aforementioned parameter ranges. The uncertainty associated with the outcome of a stochastic
32 simulation was estimated through the Monte Carlo technique and 100 simulation runs. T his number was chosen as a compromise between short computational time and high precision of confidence intervals around the mean predicted area of infestation. A set of model realizations describing the effect of chang es in parameter values on potential areas of infestation were also considered in a sensitivity analysis. Model P arameterization B asic data relevant to several termite species were used in order to parameterize the model. Unfortunately, there is not suffic ient data to calibrate the model directly against N. corniger at the Dania Beach site. The age of colonies at the first production of alates, which varies between different termite species, can be derived from related literature studies. Typically, a colon y takes four to six years from its creation to reach maturity and start the production of alates ( Collins, 1981 ) In this work the baseline value of the age of first production was set to 4 years. Lifespan es timates are approximations because they only reflect laboratory conditions. Estimated maxima ranged from 15 years old in Macrotermes bellicosus ( Smeathman ) ( Keller, 1998 ) to 20 years in Pseudacanthotermes spiniger ( Sjoestedt ) and P. militaris Hagen ( Conntable et al., 201 2 ) In this work, an age threshold of 20 years was set after which a colony dies. The maximum distance of pheromone attraction currently reported is 2.5 3 m by males ( Leuthold and Bruinsma, 1977 ) Here, the baseline value for the mode l was set at 3 meters. The density of termite colonies over a certain patch of land is related to its specific biology, ecology and behavior ( Adams and Levings, 1987 ) No specific literature sources studied the density of N. co rniger 's colonies within an urban environment.
33 However, a study found a density of approximately 7 colonies per hectare in a primary forest in Panama ( Thorne, 1982 ) which was used as a baseline value in our model. Literature sources treating the topic of alate predation or alat e flight success rate are scant. Both predation and injuries typically occur as alates start leaving the nest (i.e. pre flight), in flight (bats and birds), and as soon as they alight on the ground or on a tree (i.e. post flight) to search for a mate. Empi rical observation of alates of a different invasive termite species, Cryptotermes brevis (Walker) found an approximate survival rate of 1%, excluding predation ( Scheffrahn et al., 2001 ) Factors affecting the outcome and the success of the dispersal flight include environmental conditions, number of alates, sex ratio, proportion of alates eaten by predators, and efficiency of the post flight mating behavior ( Noirot, 1990 ; Nutting, 1969 ) A recent field study for two termite species showed that, despite the presence of 40 mature colonies over an area of one hectare producing approximately one million alates every year, no new colonies were found ( Conntable et al., 20 12 ) In this work the baseline value of the overall survival rate was set to 0.01 (1%). T his can be thought to be a realistic estimate considering all the aforementioned factors 2 Although sex ratios of termite alates are variable, they tend towards p arity (Jones et al. 1988). In N. corniger individual colonies produced alates whose sex ratio was around 1:1 ( Darlington, 1986 ; Thorne, 1983 ) Therefore, the baseline value of the prevalence of male alates in the colony was set to 0.5 (50%). Field studies ai ming to precisely assess the size of an alate crop in individual colonies are rare. Several colonies of N. corniger have been compared and a noticeable 2 Scheffrahn, personal communication
34 variation in production of alates was found. Mature colonies, whose population size ranges between 50,00 0 400,000 individuals, produced 5,000 25,000 alates ( Thorne, 1983 ) The production of alates likely depends on factors such as availability of food resources, health of queen(s), colony age, and colony specific history. All factors are not easily assessed during the short time frame given in field sampling. In another invasive termite, Coptotermes formosanus the alate production of a single colony was over 68,000. In this case, sex ratio was 1:3 (F:M) ( Su and Scheffrahn, 1987 ) related alate production was set defined as follows: (i) no production of alates until a colony reaches 4 years of age, (ii) 1,000 ala tes between 4 and 9 years of age, (iii) 10,000 alates between 10 and 14 years of age, and (iv) 100,000 alates between 15 years and the age at which a was also defined with a greater productio n of alates at an earlier age: (i) no production of alates until a colony reaches 4 years of age, (ii) 10,000 alates between 4 and 9 years of age, (iii) 50,000 alates between 10 and 14 years of age, and (iv) 100,000 alates between 15 years and the age at w hich a colony dies. This alternative scenario is tested in our sensitivity analysis (see Results and Discussion section oversimplification, it is likely to match an average magnitude that is otherwise impossible to cal ibrate with precise empirical data 3 Termite alates are weak, erratic fliers. On average, alates are capable of flying a few hundred meters on their own ( Nutting, 1969 ) Flight dista nces have not specifically been estimated for N. corniger However, it is possible to estimate this model parameter 3 Scheffrahn, personal communication
35 based on findings for other termite species. Mark recapture studies using light traps gave the first empirical measurements of termite fligh t ability A maximum distance of 892 m has been recorded for C. formosanus ( Messenger and Mullins, 2005 ) In an endemic habitat, alates may fly far enough to ensure that a mixture of different colonies is created with swarm aggregation ( Husseneder et al., 2006 ) However, for an exotic population to spread, alates fly into uncolonized areas lacking conspecifics with whi ch to mate. Recently, a new maximum distance record of about 1.3 km has been recorded by Mullins and Messenger in New Orleans, LA 4 Alates of Odontotermes formosanus ( Shiraki ) were capable of flying between 120 and 743 m ( Hu et al., 2007 ) Other studies recorded only a few dispersal flights covering about 300 m for termite species belonging to the Termitidae family ( Mill, 1983 ) to which N. corniger belongs, or 460 m for C. formosanus ( Ikehara, 1966 ) In this study, dispersal flight distances were sample d from an exponential distribution. This allows for both short and rare longer dispersal events. In a unique mark recapture study recently completed in New Orleans, LA, data collected for alates of C. formosanus irical histogram derived from several recorded flight distances (Mullins, unpublished data). T he mean of the exponential distribution was estimated based on the aforementioned empirical data and literature findings. The baseline value used as a mean disper sal distance for the simulation model was set to 200 m. Two factors that affect alate dispersal distance during the swarm season are wind velocity and light intensity. In most cases, the flight is only initiated if the wind velocity stays below 3.5 km/h ( Leong et al., 1983 ) Moreover, termites are extremely 4 Mullins, personal communication
36 prone to injuries, hence windless or low wind conditions are preferred. Given the impossibility of forecasting wind direction, wind speed, and light intensity in a multi year simulation model, alates are assumed to be able to fly in any direction and sample all angles (in radians) from a uniform distribution. M oreover, the present model was used within an urban environment, where light intensity is quite uniformly distributed and therefore it will not likely affect the model outcome. Application of M odel to S pecific S tudy A rea and D ata N. corniger was first repo rted in Florida in May 2001, in Dania Beach, Broward County, FL ( Scheffrahn et al., 2002 ) The discovery represents the first record of a non endemic and land based establishment of a higher termite (Family Termitidae) in the continental U.S. It is likely that the infestation was the result of dockside flights from an infested boa t or shipping container, probably a decade before the discovery, but no specific source was identified ( Scheffrahn et al., 2002 ) Starting in early 2003, a previously delineated area was targeted for a deliberate erad ication campaign of this invasive pest. In January 2003, an area wide visual survey was conducted for nests, foraging tubes, foraging sites, and debris harboring living N. corniger However, most N. corniger nest locations were cryptic and even an exhausti ve survey is likely to miss some infested locations, especially in the case of young colonies. In 2006, survey work was discontinued due to budget cuts before being re activated in 2011 ( Scheffrahn et al. 20 14 ). Exact sample locations were recorded using a GPS device and later imported into a database. A total area of 200 acres (approximately 81 ha) was surveyed, 20% of which had active infes tations. Several epigeal nests of different diameters were fo und at the base of both live and dead trees, in tree cavities above ground, and foraging tubes
37 extended upward of 10 m on trees and palms ( Scheffrahn et al., 2002 ) The maximum separation between active sites in north sout h and east west direction was approximately 1 km. A newly funded 2011 2012 survey revealed new infested locations. No pest reoccurrence was observed within the areas originally surveyed between 2003 and 2006 (Scheffrahn, unpublished data). Fig. 2 3 shows the known infested area in 2003, with a zoom over the recorded GPS locations of all sampled termites. The total area covers less than 0.25 km2 and consists of commercial, residential, marina, and vacant wooded properties. Results and D iscussion The stocha stic outcome of 100 computer simulations can be grouped and areas predicted as occupied by the model in at least one simulation run. Similarly, a l l areas predicted as occupied in at least half of all runs. hat are predicted as infested in all runs. Fig. 2 4 shows a snapshot of the spatial expansion of N. corniger through time as predicted by the baseline simulation model, with color coding to represent the different occupancy envelopes. Between 2003 and 2004 in the model there was a larger expansion in the areas surrounding the first surveyed locations compared to all other time frames. There are two reasons for that: (1) alates fly in all directions and therefore, if the habitat is suitable, fill in all the voids; (2) After 2004, most of the areas toward the center of invasion had already been invaded and therefore occupied by at least one colony. representing only areas that were not occupied in all simulation runs, hence they overestimate the predicted area and show a much larger extent than was likely to have
38 been invaded Areas cover ed by the "100%" envelope can be used to plan a first survey and either quarantine or eradicate the infestation. The other envelopes, instead, can be used as a "worst case scenario", thus used as a maximum perimeter to plan a more effective eradication pro gram. Overall the expansion seems to proceed slowly and it is possible to observe some barrier effect represented by both highways and the airport ground on the shape of the predicted surfaces in the East North East directions. Finally, a few isolated spot However, these spots may have been predicted by a single simulation run out of 100 and they should not be looked at as a threat. The contribution of each model parameter to the final outcome of t he computer simulation is assessed with a sensitivity analysis. This is typically done by slightly changing the value of a given model parameter while keeping the other model parameters constant. Based on the change in output one can estimate how the uncer tainty in the model output can be apportioned to uncertainty in that parameter. T he importance of each parameter was evaluated through a set of metrics, which are: covered area, absolute area growth, relative area growth. All measures are expressed as Monte Carlo (or multi run) averages, i.e., as arithmetic means of all 100 simulation runs. For six out of the seven parameters selected for the sensitivity analysis, as introduced in Table 2 1, the simulation was run with two alternative values, giving a total of 12 alternative model realizations in addition to the baseline simulation. Further, a single change of value was tested for variable SCR because only in observing the effect of a different age dependent re production structure and did not have empirical data to justify more realistic alternative scenarios on that parameter. Detailed results from the
39 global sensitivity analysis are shown in Appendix A ( Table A 1 ). Relative and absolute growth rates in the tab le refer to changes in area compared to the previous year. Here, for the sake of brevity, t he sensitivity analysis results are reported using line charts and selecting the average predicted area of infestation through time as a representative measure of ch anged parameter settings. Fig. 2 5 shows the charts for the seven tested parameters. Each chart also contains a line of the baseline model as a reference. The parameters that have the largest overall influence on the model outcome, considering all evaluati on metrics, are SCR ( scenario of amount of alates generated by a colony) SURV ( overall survival rate of alates ), PHR ( maximum pheromone attraction distance) and DIST ( mean dispersal flight distance ). The parameter MAR ( prevalence of male alates in the co lony ) has the smallest effect. Both AFP (age of first production of alates) and DEN (density of colonies /ha) have a relatively small effect. When SCR is set to "High Profile" there is a large and sudden increase in the predicted infested area after the fi rst four years, as described in the Model parameterization subs ection A higher number of alates is produced after reaching the age of first production and this increase is far more rapid compared to the "Low Profile" used in the baseline model. The PHR pa rameter has a large effect as it sets the rule for the maximum distance within which alates can find a mate. When the radial distance is reduced by two meter s the final predicted area is reduced to less than half of its corresponding baseline value. The S URV parameter controls the percentage of alates that are able to survive predation and find a mate. Therefore, the higher the percentage, the higher the chance to create new colonies at any time step. In general, the effect of a change in a model parameter ac cumulates over time. As an example, Fig. 2 6 ( B C ) shows the effect of a change in
40 the SURV parameter on the predicted area of infestation in the study area. For the sake is shown To corroborate our simula tion model, all newly infested sites that were discovered in 2012 were included Figure 2 7 (right image) shows the infested areas predicted by the baseline simulation model with all three occupancy envelopes using 2003 sample sites as seed points (left im overlaps well with the 2012 overestimate termite spread. The main goal of this work was to develop a stochastic individual based simulation model that would give regulatory agencies a tool to anticipate possible areas of infestation and, at the same time, optimize the allocation of human and financial resources toward an eradication effort. Model output may be used by either local authorities or pest control agenc ies to draw one or more areas of intervention instead of randomly inspecting an unknown perimeter with subsequent waste of resources. For example, a greater amount of economic resources could be assigned to those zones nvelope. H indcasting was used in order to predict which areas in Dania Beach, FL, would have been infested up to the present if no eradication plan had ever been implemented. The model presented in this study is a generic model for termite s and can be applied to any species after proper calibration of all the parameters. To make the model more realistic, t he complexity o f a termite invasion was captured by including several of the ecological biological characteristics that control the dynamics and speed of their natural dispersal
41 Some limitations of the model presented include the precision of the estimates used to parameterize it. In some case s parameters had to be estimated based on literature findings on termite species that are not the same as the one modeled. Unfortunately, this was necessary whenever an empirical estimate could not be found for N. corniger Although the lack of precise estimates for N. corniger may affect the final outcome of the model, all values reflect a general t endency shared by most termite species. The precision of the model presented in this study will greatly benefit from newer and better empirical estimations for the species being modeled. Whenever calibration data are missing or scant, a consultation with a termite exper t is recommended Future research may expand from our work and implement a micro level simulation model to simulate multiple dispersal steps within a single year. Moreover, future implementations may include, among other parameters, prevailing breeze direc tion and distance from city street lights for nocturnal dispersing species. The Monte Carlo technique is used to assess the uncertainty associated with the stochastic outcome of each model and to obtain an approximation of the answer to the problem O ccupa ncy envelopes were used in order to estimate areas of infestation with different likelihoods. Although the nature of the available data does not allow the use of a traditional model validation technique, th e comparison with field samples via hindcasting pr ovides at least some support to our conclusions. Results show that the areas predicted as infested in all simulation runs by our baseline model match all empirical sample locations well. A sensitivity analysis was used to check for the importance of each m odel parameter, indicating that in particular, the parameters settings for the amount of alates
42 generated by a colony, overall survival rate of alates maximum pheromone attraction distance and mean dispersal flight distance heavily influenced the final o utcome of the model. T his study is potentially beneficial to termite science, pest control agencies, and to a general audience. The simulation model was implemented using the open source R programming language. The functions are freely available to the use rs and flexible to facilitate use in different future applications The source code can be found at https://github.com/f tonini/Termite Dispersal Simulation
43 Figure 2 1. Structure of the initialization steps involved in the termite dispersal simulation model.
44 Figure 2 2 Core subprocesses involved in the individual based simu lation algorithm at any generic time step.
45 Fi gure 2 3 Location of samples of N corniger collected during a field survey in 2003 The background satellite ima ge on the top right corner was taken from a set of historical images in Google Earth.
46 Figure 2 4 Snapshot of the areas predicted as infested with N. corniger by the baseline dispersal simulation model. Yellow, orange, and red cells indicate the > 0%, >50%, and 100% occupancy envelopes, respectively. Top left map: dots represent samples of N corniger collected during a field survey in 2003, while green cells indicate the approximate areas of initial infestation.
47 Figure 2 5 Sensitivity analys is charts for the termite dispersal simulation model Each of the seven parameters is compared to the baseline simulation model (blue line). Red and green lines represent the models with a small change in a given parameter, leaving all the other variables unaltered.
48 A B C Figure 2 6 Model sensitivity to the SURV parameter. A ) Baseline simulation model. B) SURV = .005 (0.5%) C) SURV = .02 (2%).
49 Figure 2 7 E valuation of the termite dispersal simulation model Original and predicted infested areas by N. corniger with 2003 and 2012 sampled termite locations.
50 Table 2 1 Model parameters: abbreviations, definitions, and their baseline values. Parameter Definition Baseline Value Source AFP Colony age at first production of alates 4 yrs ( Collins, 1981 ) PHR Maximum pheromone attraction distance 3 m ( Leuthold and Bruinsma, 1977 ) DEN Maximum density of colonies per hectare 7 ( Thorne, 1982 ) SURV Ov erall survival rate of alates* 0.01 (1%) ( Scheffrahn et al., 2001 ) MAR Prevalence of male alates in the colony 0.5 (50%) ( Darlington, 1986 ; Thorne, 1983 ) SCR Scenario of amount of alates g enerated by a colony Low Profile (see Model parameterization subs ection) (Scheffrahn, personal communication) DIST Mean dispersal flight distance 200 m ( Scheffrahn, personal communication ) Overall percentage of alates surviving all phases of a dispersal flight
51 CHAPTER 3 STOCHASTIC SPREAD MO DELS: A COMPARISON B ETWEEN AN INDIVIDUAL BASED AND A CELL BASED MODEL FOR ASSE SSING THE EXPANSION OF INVASIVE TERMITES OV ER A LANDSCAPE The Integration of GIS with Spread Models Termites are serious pests of the urban environment and are responsible for severe damage to structural lumber ( Rust and Su, 2012 ) It has been estimated that every year, in the continental United States alone, termites cause property damage up to billions of dollars ( Edwards and Mill, 1986 ) Although human assisted (also called anthropogenic) dispersal can play a significant role in the rate of expansion of an established exotic sp ecies, its nature is complex and unpredictable. Termites are not hitchhiking insect pests, i.e. they cannot easily be transported aboard commercial vehicles such as cars or boats. More specifically, the nes ting core of a termite colony must be moved intact and both a water and food source must be available to the core throughout the movement ( Hochmair and Scheffrahn, 2010 ) The natural spread of termites is expected to proceed fairly slowly for two main reasons: (i) they are weak fliers and their reproductives (winged individuals) can fly only a few hundred meters from the parent colony on average each year ( Nutting, 1969 ) ; and (ii) a termite colony takes at least four years to mature and release the first reproductives from the nest ( Evans, 2011 ) Given the unpredictability of human assisted movements it will be more useful here to restrict our study on termite dispersal by natural means and attempt to anticipate the rate and direction of termite spread. The recent integration of Geographical Information Systems (GIS) with simulation models, coupled w ith advances in computing power, has allowed the development of spatially explicit and dynamic simulation models, such as individual based models
52 (IBMs, hereafter), and cell based models, such as cellular automata (CA, hereafter) ( Brown et a l., 2005 ) Both types of models have been introduced in ecological modeling to capture the inherent complexity of various real world problems ( Steyaert, 1993 ) as alternative approaches to solving mathematical sets of equations such as partial differential equations (PDEs) ( Alexanderian et al., 2011 ; Holmes et al., 1994 ) Moreover, the intricacy of a physical environment limits the applicability of mathematical models for modeling realistic dispersal of invasive species ( Pitt, 2008 ) GIS allows for spatial complexity and simulation modeling takes GIS visualizations into the domain of temporal dynamics. Early examples of GIS based simulation models include land use change based on the Markov process ( Burnham, 1973 ) and discrete s tate models of flow and transport (e.g. transport of pollutants) ( Maidment, 1996 ) Other more recent computer intensive approaches have been able to incorporate ecological information at a population level with landscape complexity ( Wiegand et al., 2004 ) While IBMs incorporate rules to describe the interactions between individual units, such as organisms, in which each individual can have a different set of behavioral, physiol ogical, and other properties ( Huston et al., 1988 ) in cell based modeling the basic units of the simulation are discrete spatial cells, which can transition among different states (e.g. empty, infested, quantified occupation, etc .). CAs are the most basic form of cell based model and are able to model complex dynamics from biological, social, or physical processes by using simple rules ( Toffoli and Margolus, 1987 ) Several improvements have been made to basic CAs in the past decade to include underlying landscape differences, more comp lex rules for changes of state, and stochastic parameters. Some examples of cell based models include wildfire spread
53 models ( Clarke et al., 19 94 ) spread models of invasive ants ( Pitt, 2009 ) spatial dynamics of urban and regional systems ( Clarke and Gaydos, 1998 ) and epidemic propagation ( Mikler et al., 2005 ; Morley and Chang, 2004 ) In the past years, different computer simulation approaches have been either compared ( Ajelli et al., 2010 ; Lett et al., 1999 ) or integrated ( Sudshira, 2010 ) In this work, a new stochastic cell based model was developed in order to simulate the spatiotemporal spread of invasive termites. The model developed herei n includes a set of relevant ecological parameters that can be modified according to the species of interest. A recently published stochastic IBM for the spread of invasive termites ( Tonini et al., 2013 ) was used as a benchmark for comparison with the new model in terms of both predictions and computational time. The purpose of the new m odel is to improve computational speed over the IBM model, which is costly in terms of computational time. T he model proposed herein likely represents a substantial improvement to the IBM and could give regulatory agencies a better tool for targeting speci fic areas for survey, eradication, or quarantine effort. The format of the chapter is as follows. The new model, its framework and main parameters are described in Materials and Methods Results are presented in the Results s ection, while their interpretation is discussed in the Comparison of Model Results and Discussion section. Finally, conclusions on the advantages and limitation s of the models tested herein are provided in the Advantages and Limitations of the Model sectio n. Materials and Methods A side by side comparison between the proposed cell based model and the IBM by Tonini et al. (2013) was done by establishing a common ground in order to discount
54 unwanted effects due to their different frameworks. Both models were used in two different settings: (i) a homogeneous landscape, where different sets of N termite colonies were randomly located and considered as sources of invasion, and (ii) a case study of an invasive termite, Nasutitermes corniger (Motschulsky), in Dania Beach, Florida ( Scheffrahn et al., 2002 ) The average area occupied at each time step was computed over all model replications (Monte Carlo estimation) and its value compared between the two modeling approaches. Moreover, the computational time was considered to assess the speed of the proposed model against the established IBM. Model D esign The cell based model proposed herein is spatially explicit and stochastic with the purpose of simulating the spatiotemporal spread of a termite invasion by natural means over a realistic landscape (e.g. urban environment), starting from a set of initial su rveyed colonies. The basic units of the model are cells of a grid with a pre defined spatial resolution. Each cell may assume one or more states expressed by the number of termite colonies contained in it. The temporal scale chosen for the simulations was discrete with one time step representing a year. In order to keep the model comparable to the IBM by Tonini et al. (2013), the resolution of the spatial grid was set to 100 x 100 meters (= 1 hectare). The main ecological parameters, shown in Table 3 1, we re kept identical to the IBM. A general description of each one of them follows here. Colony age at first production of reproductives The age at which a termite colony starts generating the first crop of winged reproductives (AFP, hereafter) contributes t o the rate of spread of an invasion. The earlier the first reproductives are generated, the faster the invasion will proceed.
55 Pheromone attraction distance Termite reproductives find potential mates of the opposite sex by sensing and responding to pheromo nes after concluding the dispersal flight ( Bordereau and Pasteels, 2011 ) A new colony begins with a male female (i.e. king and queen) couple of unwinged reproductives starting the nest in a proper substrate, such as soil or wood. The maximum pheromone attraction distance (PHR, hereafter) affects the chance that two heterosexual individuals find each other after the dispersal flight (see Fig. 3 1). The smaller this maximum distance is set, the smaller the number of new termite colonies will be. Colony density The maximum density of termite colonies over an area (DEN, hereafter) avoids overpopulation of grid cells in the simulation model. The higher this parameter is, the higher the chance for a nearby cell to become infested. Survival rate The overall survival rate of reproductives (SURV, hereafter) in a termite colony determines the number of reproductives that sur vive a dispersal flight and can thus potentially mate. Survival is expressed as an overall rate that considers both predation and injuries that typically occur as reproductives start leaving the nest (i.e., pre flight), volent predators in flight (capture by bats and birds), and nonvolent predators (e.g. ants and herps) as soon as they alight on the ground or on a tree (i.e., post flight) to search for a mate. Male prevalence The prevalence of male reproductives in a colony (MAR, hereafter) describes the p ercentage of males among the reproductives and affects the chance that a heterosexual pair unites within the pheromone attraction distance after the dispersal flight. That is, the chance decreases the further the value deviates from 50 percent of males.
56 Am ount of reproductives generated by a colony The number of reproductives generated by a colony (SCR, hereafter) after reaching maturity (determined by the aforementioned AFP parameter) increases with age. Different scenarios representing the total number o f reproductives by colony age can be chosen. A higher number of reproductives increases the chance of new colonies being formed. Dispersal distance The mean dispersal flight distance (DIST, hereafter) determines how far termite reproductives are able to f ly, on average, on their own. The higher this parameter value, the faster the termite spatial expansion will proceed. The age of each colony is tracked along with the number of colonies within each grid cell. Depending on the number of colonies, their age, and the SURV parameter, each cell can be expressed as the total number of reproductives available to either stay within it or fly toward neighboring cells. The main process involved in the spread of termite colonies is the dispersal of reproductives. Typi cally, dispersal flights are initiated in the beginning of the wet season and are triggered by environmental factors ( Jones et al., 1988 ; Martius, 2003 ; Nutting, 1969 ) Termi tes can only spread across a landscape and form new colonies beyond the territory occupied by the mother colony by using dispersal flight (excluding occasional and unpredictable human assisted introductions). In this study, as it was assumed for the IBM st udy, the dispersal and the subsequ ent creation of new colonies were simulated as a single annual event. This simplification is necessary because of the extreme uncertainty in the frequency and magnitude of dispersal flights during each flight season. A two parameter dispersal kernel ( Clark, 1998 ) was used to express the proportion of reproductives that fly up to a certain distance from the mother colony:
57 (3 1) where ( ) is the gamma function, is the average distance parameter, c the shape parameter, and x was set to have the same value of the baseline distance parameter in the IBM equivalent, i.e. 200 m. The kernel is exponenti al when c = 1, "exponentially bounded" when c tailed" when c c = 2. The probability kernel used in this study was chosen to be exponential to be consistent with the IBM model but can be changed to match a desired shape for the dispersal function. Probability densities were generated from a grid expressing values as Euclidean distances (in meters) from the central cell. This was accomplished first by replacing x in the equation with all the distances and then by normalizing t he resulting f(x) values with their sum in order to obtain a true probability density kernel to use as weight matrix in the focal operation. Termite reproductives find potential mates of the opposite sex by sensing and responding to pheromones after conclu ding the dispersal flight. A new colony begins with a male female (i.e. king and queen) couple of unwinged reproductives starting the nest in a proper substrate, such as soil or wood. In the IBM model, the final location of each individual reproductive was followed and mating was assumed when two individuals of the opposite sex came within a sufficiently short distance. Thus reproduction was conceptually simple in the IBM, but very costly in terms of computational time, which the cell based model attempts t o avoid. Therefore, a new approach to estimating reproduction was used in which the number of new colonies formed by N reproductives, given different values of the MAR parameter, is again
58 stochastic. This number was derived from a separate simulation study of individually moving reproductives, from which a permanent lookup table could subsequently be used for all future simulations of the cell based model. Using a lookup table rather than performing a detailed individual based simulation for each reproducti ve event saves a great deal of time in simulations. In order to account for stochasticity in the formation of new colonies, a sequence of values of N randomly placed reproductives was used and 1000 replications were run for each one of them. Fig. 3 2 shows some of the empirical histograms extracted from the aforementioned simulation for some values of N. A more extensive set of histograms can be found in Appendix B A numerical procedure was necessary because a formal mathematical derivation would be highly complex. Specifically, this is a combinatorial problem, in which any reproductive individual can be by chance a close neighbor of multiple other individuals. Moreover, once a heterosexual pair is found within a pre defined pheromone attraction radius, tho se individuals that have mated will not be available any more for others to mate with. Due to this complexity, simulations to develop a look up table for future use was the most practical solution. Because the cell based model is stochastic, it was run 100 times in order to account for the uncertainty associated with the outcome of a stochastic simulation. The number of replications was chosen to allow comparison with the IBM and was large enough to produce tight confidence intervals around the mean predict ed area of infestation.
59 Framework and M odules The simulation algorithm was implemented using a set of R functions ( R Development Core Team, 2013 ) and the source code is provided freely 5 The simulation is composed of several modules ordered in a sequential manner. The main script calls external modules that follow the steps involved in the spread of a termite invasion from the creation of reproductives to the final creation of new coloni es, which are identical to those of the IBM study by Tonini et al. (2013): (1) the ecological parameters are set; (2) in case the simulation is meant to be run over a realistic landscape, a suitability habitat GIS layer is read; (3) a set of initial coloni es (sources of invasion) are either read from an external file with spatial coordinates associated to each record or randomly located over the study area; (4) a spatial grid with pre defined resolution is created and each cell is assigned the number of col onies falling within it; and (5) if the DEN parameter was set at (1), a random elimination of colonies exceeding the limit on colony carrying capacity within each cell is performed. During step 3, in the absence of specific information, each colony is assi gned a random age between 0 and a pre defined maximum lifespan. Fig. 3 3 shows a conceptual representation of the main modules (four boxes to the right) involved in the simulation algorithm at any generic time step. At each time step, the age of each colon y is increased by one (age increase module) and colonies exceeding the maximum lifespan age are removed. After each cell has been updated with the number of colonies within it, termite reproductives are generated depending on the colonies' ages (alate gene ration module) and their total 5 https://github.com/f tonini/Termite Dispersal Simulation
60 number summed over the grid. At this stage, a dispersal kernel is used as a matrix of weights (probabilities) in combination with a summation function to calculate numbers of reproductive arriving at a focal cell from both it self and from neighboring cells. Therefore, a focal ("moving window") operation, shown in Fig. 3 4, can be used to return new values for the spatial grid representing the final number of potential reproductives available after the dispersal flight that may placed such that the focal cell is in the center. Estimates of the probabilities of reproductives from each of the eight surrounding cells and the focal cell itself arriving at the focal cell constitute the disper sal kernel. These probabilities are then used as multipliers of the number of reproductives in the focal and adjacent cells to calculate an estimate of the number of new reproductives that arrive in the focal cell. Next, the number of new colonies formed i n each cell (new colonies generation module) is determined by randomly sampling from the look up table (see Model design subs ection). If the model is run over a realistic landscape, all cells occupied by colonies falling within a patch of habitat classifie d as unsuitable are removed (habitat suitability module). Before moving to the next time step, the spatial grid is saved as a geo referenced raster grid to be overlaid at the end of the simulation with all other replications for a given year. This operatio n is the same used by Tonini et al. (2013) to compute a final occupancy envelope. Results Homogeneous L andscape Several initial scenarios were tested by randomly placing termite colonies across a spatially homogeneous study area, without any habitat limita tions for the species. These initial scenarios were used as the first time step to run both the new cell based model
61 and the IBM by Tonini et al. (2013) under the same conditions. Fig. 3 5 shows the results obtained from the stochastic outcomes of 100 mode l replications after 20 time steps for both the IBM and cell based models, considering different values of N randomly located colonies. The outcomes of the 100 replications are grouped such that they can be visualized by three color coded occupancy envelop occupancy envelope represents all areas predicted to be occupied by the model in one predicted areas that were occupied by at least half of all simulatio n runs. Finally, the all runs. Although the shapes of the envelopes computed by the two approaches are similar, there are a few differences, especially in the broadest (yell ow) envelope. Fig. 3 6 shows a comparison between both simulation approaches in terms of occupied areas over 20 years of simulation and grouped by the same aforementioned occupancy envelopes. The broadest envelope (yellow) shows the largest differences be tween the two approaches, with the cell based model predicting a larger infested area, while there C 1 (Appendix C ) shows detailed results of the multi run average occupied area at each tim e step by both modeling approaches. The computational time of the two modeling approaches needed to complete the simulation and to save all results to external geo referenced raster grids was also compared. Fig. 3 7 shows a comparison of the speed (in hours) given the different values of N randomly distributed initial colonies. All simulations were run on a 3.30 GHz Intel(R) Core(TM) i7 3960X CPU with 32 GB RAM.
62 Case S tudy: Dania Beach, Florida The proposed cell based model was also applied to the case study presented in Tonini et al. (2013) in order to compare its stochastic outcome to the IBM. Both models used surveyed N. corniger colonies ( Scheffrahn et al., 2002 ) as initial sources of dispersion and the simulated spread after 9 years was compared to the observed infested sites discovered in 2012. Values used for the Dania Beach case study were assumed identical to the ones used for the IBM. Fig. 3 8 shows a snapshot of two years extracted from the simulation in both approaches and the three types of occupancy envelopes (see Homogeneous landscape se ction) The infested locations by N. corniger observed in different years are also shown. Comparison of Model Results and Discussion Fig. 3 5 provides a comparison of the results between the cell based model and the IBM after 20 years of simulation. Although som e differences can be seen across all envelopes, the number of initial colonies does not seem to have a great impact on the differences between the spatial extents of the areas predicted by the two approaches. Visual inspection suggests that the broadest en velope ("> 0%") is the one with consistently larger differences. To assess envelope dissimilarities throughout all simulation time steps, one can refer to results in Fig. 3 6. An increasing trend in the dissimilarity between envelopes mirrors with an incre asing number of initial colonies. However, while the and the 100%" envelopes look more similar for the two modeling approaches and do not show a consistent over/under estimation over time, the "> 0%" envelope confirms the pattern observed in Fig. 3 5. Specifically, the cell based model seems to always estimate a slightly larger infested area, and this difference starts to increase approximately after the first 2 3 years of simulation. This
63 can be explained by the larger standard deviation found in the estimated occupied area over the 100 model replications for the cell based model for the "> 0%" envelope (see Table C 1, Appendix C ). Since the broadest (yellow) envelope is created to show spots that are predicted by one or more model replications, it is subject to high stochasticity, so that it is natural that a higher standard deviation in the occupied area is reflected more for this envelope than for the others. As pointed out by Tonini et al. (2013), areas predicted by both the "> 0%" and e nvelopes are conservative estimates that are likely to overestimate the true expansion and should not necessarily be looked at as a serious threat As opposed to this, areas predicted by the 100%" envelope should instead be carefully inspected and monitor ed by local authorities as they are consistently predicted to be infested by all model replications. Higher values were also tested for the DEN parameter, in the range from 5 to 10. All models had a similar outcome to that presented herein. The biggest adv antage of the new cell based model compared to the IBM can be seen clearly in Fig. 3 7, where computational speed is compared for several values of initial seed points (termite colonies). The larger the number of initial points, the longer it takes for the IBM to complete the simulation, as opposed to the cell based model which seems to tend asymptotically to a fairly stable value. Specifically, with 10 initial colonies the computational time of the IBM model already doubles that of the cell based one (~3.5 hours vs. ~1.8 hours). This difference increases further up to a factor six with 50 initial colonies (~15.6 hours vs. ~2.5 hours). This result should encourage the use of the proposed cell based model over the IBM, unless the simulation is being used to s imulate the spread of invading termites from a small set of initial surveyed colonies, in
64 which case the IBM is not much more costly and may be more accurate. Moreover, the aforementioned result was obtained from multiple simulations over a homogeneous lan dscape. Models allowing a higher maximum density of colonies per grid cell tended to slow down the simulation in both modeling approaches. However, while the cell based model only had a slight variation in speed, the IBM suffered more extreme effects, whic h almost halved its speed. The results obtained with the cell based model for the Dania Beach case study are also encouraging in the sense that both approaches have a similar outcome. Not only do all three occupancy envelopes approximately cover the same a reas, but the infested sites found between 2004 and 2012 all fall within the most probable envelope. Although this result does not represent a true model evaluation, given the limitations of the empirical data pointed out by Tonini et al. (2013), it corrob orates the cell based simulation model in so far as its results enclose the surveyed locations. The "100%" envelope is smaller in the cell based model than in the IBM because the higher standard deviation increases the chance that at least a few grid cells will not be occupied in some of the model replications. A major difference between the IBM and the new cell based model is represented by the way the underlying suitable habitat is represented, as shown in Fig. 3 9. On one hand, the IBM can use vector typ e GIS layers to consider patches of landscape that are unsuitable (e.g. roads, marshes, rivers) and all colonies that are falling within it (point in polygon operation) are removed from the simulation. On the other hand, the cell based model can only use r aster type GIS layers, hence it strongly
65 depends on the spatial resolution chosen for the simulation grid. If the grid has a coarse resolution, the underlying habitat will also be coarse and all colonies within a cell may be removed as a whole, hence losin g some important detail necessary at a small spatial scale. Advantages and Limitations of the Model The new cell based model presented in this study is a generic model for termites and can be applied to any species upon modification of the values of its ec ological parameters. Although individual based models are able to provide a very detailed scenario and incorporate both individual differences and local interactions among organisms, their computational cost is often higher compared to cell based models ( Ajell i et al., 2010 ) Therefore, whenever the spatial scale is large (e.g. state, country level) or the temporal window is long, it is recommended to use a different simulation approach. Given the findings of this study, the proposed cell based approach not only gives similar results to the IBM equivalent, but has a significant advantage in terms of computational speed. This has been shown to be particularly true when there is a large numb er of initial termite colonies at the start of the simulation. However, given the higher variability found in the results of the proposed model, it is recommended to use more replications whenever possible to have a greater precision in the perimeters deli neated by the occupancy envelopes. The degree of precision of the estimates used to parameterize the model represents a limitation of the proposed approach. As already pointed out by Tonini et al. (2013) for the IBM model, the cell based model will also gr eatly benefit from newer and better empirical estimations for a given species. It is often hard to find a precise
66 estimate for cryptic species like termites and a consultation with a termite expert is always encouraged to calibrate according to the desired species. N ewer and better empirical estimations of the ecological parameters for a given species will be of great benefit to the simulation model. Finally, it is important to remark that the least probable envelopes should only be used for delineating a g eneral perimeter of intervention, close to a "worst case" scenario, instead of randomly inspecting unknown areas. A greater monitoring, quarantine, or eradication effort could be used in those areas encompassed by the most probable envelope.
67 Fig ure 3 1 Schematic example of termite reproductives and the pheromone attraction radius within which a heterosexual individual must be found to mate after a dispersal flight. The male at the top is beyond detection range of a female while the males at the bottom left are in female detection range
68 Figure 3 2 A sampling of empirical histograms for the number of new termite colonies formed by N reproductives (alates). Histograms are derived from 1000 replications.
69 Figure 3 3 Conceptual diagram of the main modules used in the simulation algorithm at any generic time step
70 Figure 3 4 Schematic representation of the focal ("moving window") operation. Focal values are calculated for the neighborhood (dashed line) of focal cells (red in this example) using the dispersal kernel as matrix of weights, in co mbination with a summation function.
71 A. B. C. Figure 3 5 Simulation outcome comparison between cell based and individual based models A) I nitial scenarios of N randomly located colonies. B) Final outcome of the cell based model ; and (C ) final outcome of the individual based model.
72 Figure 3 6 Predicted areas (in hectares) of infestation for the cell based model (solid line) and the individual based model (dashed line) for different values of N randomly distributed initial colonies. Occupancy envelopes are color coded: ), red).
73 Figure 3 7 Computational time (in hours) of the cell based model (solid line) and the individual based model (dashed line) for different values of N randomly distributed initial colonies.
74 A. B. C. D. Figure 3 8 Snapshot of the areas predicted as infested with N. corniger by the simulation models with sampled termite locations in 2003 and the period 2004 2012. ( A B ) Individual based model. ( C D ) Cell based model. Yellow, orange, and red cells indicate the > occupancy envelopes, respectively.
75 A. B. Figure 3 9 Habitat representation in the IBM (left) and the cell based model (right). A) Termite colonies (dots) falling within patches of unsuitable habitat (in grey) are removed. B) Cells occupied by one or more termite colonies (black) overlaying cells of unsuitab le habitat (grey) are removed.
76 Table 3 1 Model ecological parameters used in both modeling approaches and their values. Parameter Definition Values AFP Colony age at first production of reproductives 4 yr PHR Maximum pheromone attraction distance 3 m DEN Maximum density of colonies per grid cell 1 (homogeneous landscape setting) 7 (Dania Beach case study) SURV Overall survival rate of reproductives 0.01 (1%) MAR Prevalence of male reproductives in the colony 0.5 (50%) SCR Scenario of amount of reproductives generated by a colony DIST Mean dispersal flight distance 200 m Overall percentage of alates surviving all phases of a dispersal flight
77 CHAPTER 4 PREDICTING THE GEOGR APHICAL DISTRIBUTION OF TWO INVASIVE TERM ITE SPECIES FROM OCCURRE NCE DATA Species Distribution Models growing interest during the past two decad es given the importance of monitoring endangered or invasive species and understand a species' response to different environmental conditions ( Guisan and Thuiller, 2005 ) Such models are often referred to as habitat models, ecological niche models, or species distribution models (SDMs) ( Elith and Leathwick, 2007 ) and h ave been applied to a variety of fields such as ecology, conservation, and biogeography. SDMs attempt to model the species environment relationships by using sites of known occurrence (presence data) and, sometimes, non occurrence (absence data) together w ith environmental variables recorded over the whole study area. In most cases, records from atlases, herbaria, or museum databases only contain information on a species' incidental observations ( Franklin, 2009 ) A fundamental limitation of presence only datasets is that the prevalence of a species, i.e. the proportion of occupied sites across the study a rea, is unknown. In recent years, several statistical methods have been proposed for modeling this type of datasets, such as inhomogeneous Poisson process (IPP), ( Chakraborty et al., 2011 ; Warton and Sheperd, 2010 ) and maximum entropy (Maxent) ( Phillips et al., 2006 ; Phillips et al., 2004 ) Other approaches use presence absence models by treating random samples chosen from th e region of interest (background samples) as absences (also called "pseudo absences") ( Elith et al., 2006 ) However, this assumption has been shown to have substantial errors ( Warton and Sheperd, 2010 )
78 In this work, a recently developed Bayesian logistic regression model adjusted for presence only data ( Divino et al., 2013 ; Divino et al., 2011 ) and the widely used maximum entropy approach were used to predict the current and future biotic distributions of two major invasive termite pests over Florida: the Asian subterranean termite (AST), Coptotermes gestroi (Wasmann), and the Formosan subterranean termite (FST), Coptotermes formosanus Shiraki. The Bayesian approach used herein has only bee n tested on artificial data prior to this study. Both AST and FST are considered among the most destructive termite species, causing severe damage to structural lumber ( Rust and Su, 2012 ) AST is endemic to southeast Asia and it is currently found mostly in tropical areas ( Li et al., 2009 ) FST is probably endemic to southern China and is found prima rily in subtropical and temperate regions ( Li et al., 2009 ) AST and FST occur sympatrically only in Taiwan, Florida, and Hawaii ( Li et al., 2010 ) AST was first found in 1996 (Dade County) and is a more recent invasive species compared to FST, discovered in 1980 (Broward County) ( Scheffrahn, 2013 ) Both species are now well establishe d pests in Florida. Regional predictions of the potential habitat of the two termite species under both current and future climate scenarios are currently lacking in the available literature. A single recent study attempted to predict the ecological niche of AST on a global scale using mostly coarse precision occurrence data derived from the literature ( Li et al., 2013 ) However, the reliability of such predictions was strongly affected by the excessive extent of the study area used for both model calibration and estimation, given the small amount of available occurrence data.
79 The format of the chapter is as follows. The study area, data, variables, and modeling approaches used are described in Materials and Methods Results and their interpretation are presented in the Results se ction and a final discussion on the advantages and limitations of the models tested herein is provided in the Discussion s ection. Materials and Methods Study A rea and S pecies D ata Florida was selected as a common study area for both AST and FST in order to compare the performance of two different statistical approaches under the same environmental conditions. Termite collection localities were taken from the University of Florida Te rmite Collection at the Fort Lauderdale Research and Education Center. Geographical coordinates of 280 and 411 separate land based colonies of AST (1996 2012) and FST (1985 2012), respectively, were used in this study (Fig. 4 1). Boat infestations ( Scheffrahn and Crowe, 2011 ) were excl uded. All database samples were predominantly collected from or near human habitations by R.H.S., pest control professionals, property owners, entomologists, and others interested in species level identification. The study area was divided into roughly 38, 000 2 km grid cells and all termite observations falling within a given cell were aggregated to a single point. In this work, grid cells were considered as independent and the probability of presence was modeled for each one of them. The chosen spatial res olution attenuates some of the bias caused by spatial dependence between nearby occurrences because termite reproductives from a mature colony fly only a few hundred meters during their annual dispersal flights ( Nutting, 1969 ) Moreover, it is a common approach to consider a sampling unit whose
80 size is equal (or close) to the resolution of available environmental data ( Elith and Leathwick, 2009 ) Finally, a 2 km resolution grid makes sure that a sufficient number of observations is available for use with statistical models. After aggregation, a total of 65 and 160 occurrences were considered for AST and FST, respectively. Predictor V ariables A set of gridded climate and environmental variables was selected (Table 4 1) based upon its ability to directly influence the ecophysiology of both AST and FST ( Gautam and Henders on, 2011 ) Data for historical climatic conditions were extracted from two sources: (i) the PRISM Climate Group database ( Daly et al., 2002 ) and (ii) t he WorldClim (1950 2000) database ( Hijmans et al., 2005 ) General annual trends such as annual total precipitation (prec), aver age daily mean dew point temperature (dew), and average daily maximum (tmax) and minimum (tmin) temperatures were obtained from the PRISM database, representative of average historical conditions for the years of available occurrence records of both AST an d FST. Two bioclimatic variables representing extreme or limiting factors such as maximum temperature of the warmest month (bio5) and minimum temperature of the coldest month (bio6) were chosen from the WorldClim database, representative of 1950 2000 avera ge historical conditions. Both WorldClim and PRISM data were obtained at a 2.5 arcmin (~4 km) resolution and further resampled down using bilinear interpolation to maintain the higher data resolution of the reference spatial grid over the study area. The a vailable time series of historical climate data from PRISM (1895 Present) allowed us to extract those years that matched historical occurrence records for both AST and FST exactly. In addition to climate variables, the U.S. Geological Survey National Land Cover Database (NLCD) 2006 (http://www.mrlc.gov/nlcd2006.php) was also considered, which
81 has a native resolution of 30 m. The database comes with 20 land cover classes, which were modified according to the following steps: (1) reduction from 20 to 8 main l and cover classes according to the NLCD 2006 product legend; (2) creation of single layers for each land cover class from the previous step; and (3) aggregation of each land cover layer from 30 m to our 2 km reference grid by expressing each cell value as the percentage of land cover contained within. Finally, centroids of grid cells occupied by termite locations were also used in some of the statistical models in order to account for the geographic proximity between collection sites across the geographic space. Locations were expressed by their projected easting and northing values. All layers, including the 2 km reference grid, were mapped using the NAD83 / Florida GDL Albers projection to minimize distance distortions throughout the study area. Most pred ictor variables in our dataset were highly correlated and their simultaneous presence in statistical models has been proven to cause several problems (multicollinearity, Farrar and Glauber 1967 ) Therefore, an a priori choice of variables was carried out in order to exclude pairs of highly correlated variables (r 0.8). A different set of models was also estimated using principal components obtained from the full set of predictor varia bles considered herein. Principal component loadings are shown in Appendix D, Table D 1 Future S cenarios The Intergovernmental Panel on Climate Change (IPCC) Fourth Assessment Report (AR4) describes a set of alternative CO2 emissions scenarios grouped under four main narrative storylines (http://www.ipcc.ch/). In this study, the A2 emission scenario was focus ed on, which is at the higher end of all the future scenarios. The
82 average increase in global surface temperature by 2100 predicted by a sophisticated set of Atmospheric Oceanic Global Circulation Models (GCMs hereafter) amounts to about 3.4C ( ranging fro m 2.0C to 5.4C) This scenario was preferred over others in order to assess the impact of a larger climate change on both termite species' potential distributions and consider it as a benchmark "worst case" scenario. Given the uncertainty associated with the path of future climate change, average projections of annual precipitation and minimum/maximum temperatures for the years 2040 2069 (referred to as 2050s hereafter) were extracted from the Climate Change, Agriculture and Food Security (CCAFS) web port al (http://www.ccafs climate.org/data/). Projections of annual mean dew point temperatures were not available from any data provider, hence this predictor variable could not be considered for future scenarios. The three following GCMs ( Diniz Filho et al., 2009 ) statistically downscaled usi ng the so called delta method ( Ramirez Villegas and Jarvis, 2010 ) were selected for th e A2 emission scenario and the 2050s time frame: GFDL CM 2.1, NCAR CCSM 3.0, UKMO HadCM3. A projected population growth scenario in 2060 was obtained from the University of Florida GeoPlan Center (http://www.fgdl.org). The dataset assumes no further popul ation growth in areas currently urbanized. First, a change in the climatological variables was assumed under the A2 emission scenario given by the three selected GCMs, assuming no change in population. Then, a population growth scenario was added together with climate change, resulting in a total of six future scenarios for each species. To create
83 "consensus" maps of projected probabilities in the 2050s, predictions were averaged over the three GCMs. Modeling A pproaches In this study, several models were c onsidered using two statistical approaches for presence only data: (i) maximum entropy ( Phillips et al., 2006 ) ; and (ii) a Bayesian linear logistic regression adjusted for presence only data, named Bayesian for presence only data (BPOD) hereafter ( Divino et al., 2013 ; Divino et al., 2011 ) The former, was presented in Phillips et al. (2004) and it is widely used for modeling distributions of species ( Elith et al. 2011 ) The BPOD, following the work of Ward et al. (2009), considers a linear logistic regression adjusted for presence only data in a Bayesian framework. Although different in their theoretical backgrounds, both methods use the Bayes' rule as an impo rtant point to calculate the probability of presence of the species conditioned on the environment. An outline of the theory, main assumptions, and modeling settings used in both approaches follows. Maximum e ntropy a pproach Maximum entropy (Maxent hereafte r) is a machine learning method that uses species occurrences and a random sample of background environmental data over a region of interest to predict species distributions. Let us define to be the probability distribution of covariates, i .e. environmental variables, across locations where the species is observed ( Y = 1), and to be the probability distribution of covariates where the species is absent ( Y = 0). The quantity of interest is the probability of presence of a spec ies, conditioned on a set of
84 environmental covariates X Maxent considers the modeling of and uses the Bayes' rule to estimate the sought conditional probability distribution: (4 1) The core of the Maxent "raw" model output is the estimate of the ratio This is accomplished by seeking an estimate of that is consistent with available occurrence data. Among several possible distributions, one that maximizes the entropy of or, in other words, minimizes the relative entropy of with respect to (measured using the Kullback Leibler divergence) is chosen. The distribution of maximum entropy, i.e. closest to the uniform probability distribution or most spread out, is estimated while being subject to a set of constraints imposed by the information available from the environmental conditions where the species occurs. Environm ental variables or functions thereof are known as "features" and are treated as an expanded set of variables to be added as terms in the model specification. A random sample of background locations informs the model about The set of constraints o n ensures that empirical averages of each feature approximate their averages at sites where the species is present (or a random sample thereof). The probability distribution of maximum entropy is a Gibbs distribution, which has an exponential form ( Della Pietra et al., 1997 ) Raw exponential values estimated by the model are scale dependent, e.g. they can be extremely s mall if the study area is large, and only represent a measure of relative suitability of each site. However, the model can also be transformed from an exponential family model into a logistic model, thus
85 making it more comparable with other machine learnin g or generalized linear/additive models ( Phillips and Dudik, 2008 ) To calculate the final conditional probability of o ccurrence knowledge of the prevalence of the species i.e. the proportion of occupied sites across the study area, is required. However, is unknown with presence only data ( Ward et al., 2009 ) In this case, the maximum entropy approach sets this quantity arbitrarily to 0.5. Bayesian for p resence only d ata (BPOD) a pproach The approach introduced in Divino et al. (2011) and fully formulated in Divino et al. (2013) cons iders a stratified case control design adjusted for presence only data developed in a Bayesian setting. When dealing with presence only data, sampling from the reference population of locations cannot be performed under the traditional random sampling desi gn. Specifically, while a random sample of presences is available, a random sample of absences cannot be obtained. Therefore, a random sample of "contaminated controls", i.e. a random sample of locations from the whole reference population (background samp le) that can also include some occurrences of the species, is matched with the aforementioned random sample from the available occurrence data ( Lancaster and Imbens, 1996 ) Building upon the approach presented in Ward et al. (2009), the B POD implements a linear logistic regression model adjusted for presence only data within a Bayesian framework. In order to estimate the regression parameters, a two level scheme is used: (1) a first level describing the probability law that generates the p opulation data; and (2) a second level using a stratified case control design, modified
86 for presence only data to select samples from the population. In a traditional logistic regression, the response variable Y = 0 marks the absence of an attribute of int erest in the population, while Y = 1 marks the presence of the same attribute. The key point in the BPOD approach is the introduction of a stratum variable Z considered as the only observable variable. Specifically, Z = 0 means that a location is collecte d from the whole reference population, while Z = 1 indicates that a location is collected from the sub population of presences. Z = 1 implies that Y = 1, while Z = 0 implies that Y is an unknown value The introduction of the stratum variable Z al lows us to define a linear logistic regression, adjusted for presence only data. Denoting by the probability that a location is sampled ( C = 1) from the set of locations where the species of interest is present ( Z = 1) and with covariates X = x the linear logistic model for presence only data can be defined as: (4 1) where q is a correction term, depending on the number of presences truly observed and the unknown number of presen ces hidden in the sample of "contaminated" controls. An approximation of q can be derived iteratively within the estimation algorithm. After prior distributions are defined over the parameters of interest (the linear coefficients and the unobserved resp onses in the sample of "contaminated" controls), Bayesian inference can be carried out through Markov Chain Monte Carlo (MCMC) techniques ( Robert and Casella, 2004 ) In particular, an algorithm including a data augmentation step ( Tanner and Wong, 1987 ) is used to obtain an estimate of the of the logistic model.
87 Evaluation of m odel p erformance In this study, model performance is evaluated according to three criteria: (i) prediction accuracy of occurrence data, i.e. model sensitivity expressed by the percentage of correctly predicted occurrences in the sample; (ii) goodness of fit, using both the Akaike Information Criterion (AIC, Akaike, 1974 ) and its corrected version (AICc, Burnham and Anderson, 2002 ) ; and (iii) ecological realism, i.e. assessing predictions against prior biological knowledge of a species. AIC and AICc for all Maxent models were calculated u sing the ENMTools ( Warren and Seifert, 2011 ) which uses Maxent "raw" suitability scores, i.e. exponential values standardized over the study area. S everal other traditional statistical evaluation metrics such as Cohen's Kappa ( Cohen, 1960 ) or the area under the receiver operating characteristic curve (AUC, Hanley and McNeil, 1982 ) are commonly used with pres ence absence (or pseudo absence) data. However, in this study no assumption s of pseudo absence for background data are made While model sensitivity was compared across all models and both statistical approaches, AIC and AICc values were only used to compa re the relative quality of each model within the same statistical approach in order to provide a from BPOD, hence values of both AIC and AICc cannot be compared across models considered in both approaches. Sampling s cheme The following background sampling schemes were considered with respect to Maxent and BPOD modeling approaches. Each Maxent model was run 16 times, with the background sample size set to 10,000 randomly selected points. Although there are not set guidelines regarding the
88 ideal number of background points to use in each situation, from a statistical point of view it is desirabl e to have a sample that well represents the domain of the model variables over the study area. Moreover, some studies found that predictive accuracy of Maxent was best with about 10,000 points (Barbet Massin et al. 2012). All other settings in the MaxEnt s oftware have been used with their default values. Each BPOD model was run 500 times, with sample size set according to the presence/background ratio of 1:4. Specifically, in AST a sample of 65 observed presences was matched with a background sample of 654 =260 locations (total sample size n=325), while for FST a sample of 160 observed presences was matched with a background sample of 1604=640 locations (total sample size n=800). The MCMC algorithm with data augmentation used 15,000 iterations (10,000 burn in) to estimate the unknown model parameters. The reason for using different sampling schemes between Maxent and BPOD is due to the fact that the two approaches have different requirements for reaching robust parameter estimates. Specifically, Maxent need s a large background sample, while BPOD needs a large number of model replications. Given these constraints, model settings were chosen E 1, Appendix E, and Table F 1, Append ix F) In both approaches, parameter estimates were obtained as averages over all model replications. Results Several models were run to predict the current potential distribution of both AST and FST. A list of the best performing models is shown in Table 4 2 for AST and Table 4 3 for FST, together with their evaluation metrics.
89 For AST, the model that reached the highest overall performance in the maximum entropy approach was M1, while in the BPOD approach it was MPC3, which used the first three principal components as covariates. The BPOD approach outperformed Maxent consistently in all models, reaching a higher sampling sensitivity. The highest sensitivity reached by any Maxent model is ~61%, hence much lower compared to the best BPOD model (~97%). Fig. 4 2 ( A B ) shows the current potential distributions of AST predicted by the best overall models in both approaches, thus BPOD MPC3 and Maxent M1, respectively. Southeastern Florida and the Keys Islands show a much higher suitability compared to other areas, matching the general pattern of recorded occurrences. Low probabilities are also predicted along the east coast up to central Florida and on the west coast around urban areas such as Ft. Myers and Tampa. For FST, the model that reached the highest overall performance in the maximum entropy approach was MPC3, while in the BPOD approach it was MPC6, using the first three principal components and all six principal components as covariates, respectively. As for AST, the BPOD approach performed consistently bet ter than maximum entropy across all models. The highest sensitivity reached by any Maxent model is ~66%, compared to ~74% for the best BPOD model. Fig. 4 3 ( A B ) shows the current potential distributions of FST predicted by the best overall models in both approaches, thus BPOD MPC6 and Maxent MPC3, respectively. Highest suitability values are associated with urbanized areas across the entire state. Although no oc currences were recorded in some urban areas, a medium to high suitability is predicted for the species in areas such as North West Florida around Pensacola, along
90 the west coast in Sarasota and Port Charlotte, along the east coast in Melbourne and Palm Coa st, and all the Keys islands south west of Key Largo. Future predicted probabilities of presence were derived using a model from the BPOD approach for both AST and FST. Due to data availability (see Future scenarios subsection ), the BPOD M1 model was chose n to predict their future distributions. Fig. 4 4 ( A D ) shows the contemporary predictions calculated using model BPOD M1 for AST and FST, respectively. Predictions are not much different from the best models that used principal components as covariates, with the exception of a few areas for FST such as the Keys Islands or the west coast of Florida where the suitability is slightly lower. Fig. 4 4 ( B E ) shows average "consensus" projected probabilities, i.e. averaged over the three GCMs, for AST and FST, respectively, under climate change conditions for the 2050s time period and given no change in land cover. The climate variables that are projected to the future from M1 are precipitation, bio5, and bio6. Urban areas in southeast Florida seem to have an in creased predicted probability of presence for AST, while for FST changes in suitability are less noticeable. The population growth scenario (Fig. 4 4 C F ) increases the percentage of areal units occupied by developed areas, thus increasing the variable "d evel" (refer to Table 4 1) in our model. The effect of a combined change in climate and developed areas increases the predicted probabilities of presence for both AST and FST. However, the effect is much more noticeable for the latter across the whole stud y area. Discussion The performance of the BPOD approach on both species was shown to be consistently better than the widely used maximum entropy approach in terms of sampling sensitivity, with few exceptions (see Table 4 3). Moreover, the best BPOD
91 model g ave more realistic predictions from an ecological perspective compared to the best Maxent model for both species. Specifically, for FST maximum entropy tends to over predict areas across the entire state, far apart from recorded occurrences, and under pred ict areas close to them (Fig. 4 3). Although this phenomenon is less pronounced for AST, areas in the metropolitan southeast Florida are under predicted nearby recorded occurrences. The BPOD approach seems to make a better use of the information from PCA derived variables compared to maximum entropy, as its predictive power increases until reaching an optimum in terms of sensitivity and both information criteria (see Tables 4 2, 4 3). However, such models behave in a slightly different manner between the t wo species. In particular, BPOD models for AST reach an optimum with a smaller number of PCA derived variables than FST (3 vs. 6 principal components, respectively). This probably means that the original environmental variables enclosed in the first three principal components are sufficient to explain the ecological niche of AST in Florida. FST is a more generic species, i.e. spread across nearly the entire study area and exploring a wider range of values in the environmental variables chosen for the study. This result also suggests that some environmental factors influencing the habitat of FST may be missing. Maxent models reported in this chapter were estimated by fitting linear responses to relationships between response and predictor variables in order to keep comparability between the two different statistical approaches. Maxent models fitting more complex responses were also tested but had a much lower predictive performance compared to the ones fitting linear features. A major advantage of the BPOD ap proach over
92 maximum entropy is that the MCMC algorithm does not require the a priori knowledge of the population prevalence as it is considered as a further parameter in the model. This overcomes the issue of prevalence specification pointed out by Ward et al. (2009). A Bayesian modeling framework allows flexibility in the treatment of uncertainty while making full inference on the probability of presence possible. However, a more formal comparison between the two statistical approaches based on artificial data is suggested for future studies. In this work statistically downscaled climate projections for the 2050s were preferred over dynamically downscaled projections, such as the CLARReS10 dataset for the Southeast United States (http:// tiny.cc/0gjd5w ). Although the latter are able to incorporate regional scale processes their spatial resolution (~10 km) was too coarse to assess the effect of variation in climate and urbanization on the same scale used for contemporary predictions for both termite species. Climate change under the A2 scenario for the 2050s has a moderate effect on both species' geographical distribution. Conversely, a combined effect of climate change with a population growth scenario has a larger impact on their p rojected probabilities, especially for FST. This suggests that both termite species are influenced by urban development much more than by climate alone. Two issues not addressed in this work are the residual autocorrelation that may still persist among nei ghboring occurrences and the problem of observer bias ( Syfert et al., 2013 ) While an attempt to reduce the first issue was done both by using a spatial resolution at which termite occurrences could be assumed independent of each other (see Study a rea and s pecies d ata subs ection) and by using spatial coordinates as
93 predictor varia bles, true spatial models may be able to improve the final predictions. However, this goes beyond the purpose of this study. The second issue would be extremely complicated to specify in the models developed herein, as data comes from different sources and involves multiple data collectors ( see Study a rea and s pecies d ata subs ection ). Moreover, the data used in this study were not collected based on road accessibility. The treatment of such a complex issue is deferred to future work.
94 Figure 4 1 Florida occurrences of AST (green) and FST (purple).
95 A. B. Figure 4 2 Current predicted prob abilities of presence for AST. A) BPOD MPC3 model. B ) Maxent M1 model. Darker red areas correspond to areas with higher probabilities.
96 A. B. Figure 4 3 Current predicted prob abilities of presence for FST. A) BPOD MPC6 model. B ) Maxent MPC3 model. Darker red areas correspond to areas with higher probabilities.
97 A. B. C. D. E. F. Figure 4 4 Current and average projected probabilities of presence for the 2050s time period. A ) BPOD M1 contemporary predictions for AST. B ) BPOD M1 projected predictions for AST averaged over the GFDL CM 2.1, NCAR CCSM 3.0, and UKMO HadCM3 global circulation models under the A2 emission scenari o. C ) BPOD M1 projected predictions for AST averaged over the GFDL CM 2.1, NCAR CCSM 3.0, and UKMO HadCM3 global circulation models under the A2 emission scenario and population growth. D ) BPOD M1 contemporary predictions for FST. E ) BPOD M1 projected pred ictions for FST averaged over the GFDL CM 2.1, NCAR CCSM 3.0, and UKMO HadCM3 global circulation models under the A2 emission scenario. F ) BPOD M1 projected predictions for FST averaged over the GFDL CM 2.1, NCAR CCSM 3.0, and UKMO HadCM3 global circulatio n models under the A2 emission scenario and population growth. Darker red areas correspond to areas with higher probabilities.
98 Table 4 1 Environmental variables used and their data source for both AST and FST. Database Variable Description PRISM prec Annual total precipitation dew Average daily mean dew point temperature tmax Average daily maximum temperature tmin Average daily minimum temperature WORLDCLIM bio5 Maximum temperature of the warmest month bio6 Minimum temperature of the coldest month NLCD 2006 water Open water or permanent ice/snow cover devel concrete, buildings, etc.). barren Bare rock, gravel, sand, silt, clay, or other earthen material, with little or no "green" vegetation forest Tree cover (> 6 m tall). Tree canopy accounts for 25% to 100% of the cover shrub Natural/semi natural woody vegetation with aerial stems (< 6 m tall) herb Natural/semi natural herbaceous vegetation (75% 100% of the cover) cultiv Herbaceous vegetation that has been planted or is intensively managed for the production of food, feed, or fiber (75% 100% of the cover) wetlands Soil or substrate is periodically saturated with or covered with water
99 Table 4 2 List of main models used for AST, their sampling sensitivity, and information criteria. M1: X (easting), prec, bio5, bio6, all land cover variables. M2: prec, bio5, bio6, all land cover variables. M3: X (easting), Y (northing), prec, bio5, all land cover v ariables. MPC x : x stands for the number of principal components used as covariates. The best models are highlighted in bold. Approach Model Sampling Sensitivity AIC AICc Maxent M1 0.61 823.4 825.2 M2 0.61 828.7 830.6 M3 0.60 857.6 860 MPC1 0.40 1066 1066 MPC2 0.61 942 942 MPC3 0.60 828.8 829.2 MPC4 0.60 829.1 829.8 MPC5 0.58 831.7 832.5 MPC6 0.57 830.2 831.1 BPOD M1 0.96 62.2 38.2 M2 0.96 63.2 38.4 M3 0.96 67.7 41.1 MPC1 0.76 118.3 112.4 MPC2 0.90 78 70.2 MPC3 0.97 48.8 39 MPC4 0.97 50.1 38.3 MPC5 0.97 51.4 37.8 MPC6 0.97 52.5 36.9
100 Table 4 3 List of main models used for FST, their sampling sensitivity, and information criteria. M1: X (easting), prec, bio5, bio6, all land cover variables. M2: prec bio5, bio6, all land cover variables. M3: X (easting), Y (northing), prec, bio5, all land cover variables. MPC x : x stands for the number of principal components used as covariates. The best models are highlighted in bold. Approach Model Sampling Sensitiv ity AIC AICc Maxent M1 0.55 2712.9 2714.2 M2 0.55 2713.1 2714 M3 0.57 2712 2713.2 MPC1 0.66 1177 1177.1 MPC2 0.59 1047.4 1047.6 MPC3 0.66 956.8 957.2 MPC4 0.57 968.5 969.1 MPC5 0.61 963.6 964.6 MPC6 0.58 961.1 962.5 BPOD M1 0.73 391.1 366 M2 0.73 391.4 367 M3 0.73 391.5 367.2 MPC1 0.05 689.2 683.2 MPC2 0.53 529.7 521.8 MPC3 0.71 396.3 386.4 MPC4 0.71 395.1 383.2 MPC5 0.72 388.6 374.8 MPC6 0.74 382.8 367
101 CHAPTER 5 CONCLUSIONS The two simulation models developed in this study and presented in Chapter 2 and Chapter 3 showed different approaches that can be used to replicate the geographical spread of invasive termites over a landscape. Although the IBM is able to capture individu al differences and represents a more refined approach to simulation over smaller spatial extents, the cell based model gave similar results and had a significant advantage in terms of computational speed. Given the stochastic nature of both simulation appr oaches and the subsequent variability in the model outcomes, it is recommended to use a sufficiently large number of replications The application of both the IBM and the cell based model to the case study of N. corniger in Dania Beach, FL, showed that the areas predicted as infested in all simulation runs of a baseline model cover ed the spatial extent of all locations recently discovered Although the nature of the available data did not allow the use of traditional validation techniques, the comparison of model outcomes with field samples via hindcasting provided some support to the goodness of both simulation approaches. The precision of the estimates used to parameterize both simulation models is crucial in order to obtain reliable prediction, hence it is strongly recommended to consult with a termite expert should any empirical value be missing in the related literature. Moreover, given the scant information available on termites, as opposed to other social insects such as ants or bees, it was in some case necessary to use parameter values from termite species different from the one modeled. A sensitivity analysis run on the IBM gave some insight on the importance of each model parameter. Specifically, the outcome of the simulation was particularly sensitive to the amount of alates generated
102 by a colony, the overall survival rate of alates, the maximum pheromone attraction distance, and th e mean dispersal fl ight distance Both simulation models were developed with the goal of keeping them as generic as possible and applicable to any termite species after proper calibration of all the ecological parameters involved. If the spatial extent of the study area over which the simulation needs to be run is large (e.g. state, country) or the temporal window is long (e.g. more than 15 20 years) the use of the cell based approach is strongly encouraged Instead, it is advisable to use the IBM model if a more detailed local simulation is preferred. The models presented herein represent a first attempt to predict and delineate areas of infestation by invasive termites from a set of surveyed locations (or from one or more points of introduction) The st ochastic component allows a more flexible and realistic representation of the dynamic of a spread over a landscape while giving different likelihoods of infestation in output. A greater amount of economical and human resources should be assigned by state a nd local regulatory agencies to those areas with the highest likelihood of infestation. The least probable areas, instead, should only be looked at as a general perimeter of intervention in order to avoid a waste of resources that would otherwise be random ly assigned in order to survey unknown locations. The statistical models described in Chapter 4 attempted to predict the current and future biotic distributions of two major termite pests, C. gestroi and C. formosanus ov er Florida. Not only the BPOD appr oach performed consistently better compared to maximum entropy models but it also made a better use of the information from PCA derived ecological variables in both termite species. Moreover, the use of a Bayesian
103 framework overcome s the issue of prevalenc e specification given by presence only da tasets of a species' occurrence by consider ing it as a further parameter in the model. However, a more formal comparison between the two statistical approaches is deferred to future work. The Bayesian framework used in this study can be applied to any termite species, but more generally to any datasets where only occurrences of a species are known. Future scenarios considered both climate change and population growth projections The effect of climate change under th e pessimistic A2 emission scenario and three different statistically downscaled GCMs for the 2050s time frame had a much more m oderate effect on both species' geographical distribution compared to the combined effect of climate change and population growth This result suggests that both C. gestroi and C. formosanus are affected by urban development much more than by climate alone. Future work should focus on modeling the residual autocorrelation among neighboring species occurrences and the problem of obse rver bias ( Syfert et al., 2013 ) The type of statistical models presented in this study can be used by g overnment regulators and the pest control industry to anticipate the arrival of invasive termite pests or predict their potential habitat under either current or future conditions
104 APPENDIX A SENSITIVITY ANALYSIS Table A 1 Sensitivity analysis results extracted for some years of the simulation. All values in the table are computed as differences between a given model and the baseline model. Positive numbers indicate that a given model has a higher value compared to the baseline. As an example, an absolute growth of 0.09 for year 2006 (AFP=2) means that the increase in the area between 2005 and 2006 for the AFP=2 model is 0.09 km 2 / yr larger than for the baseline model (AFP=4). Parameter Variable* 2004 2006 2008 2010 201 2 AFP=2 Area (km 2 ) 0.01 0.1 0.19 0.19 0.22 Abs growth (km 2 / yr) 0.01 0.09 0.03 0.01 0.02 Rel growth (%) 3 6 1 0 0 AFP=6 Area (km2) 0.02 0.05 0.09 0.14 0.12 Abs growth (km2 / yr) 0.02 0.01 0.03 0.01 0.01 Rel growth (%) 7 0 2 1 0 PHR=1 Area (km 2 ) 0.61 0.72 0.73 0.8 0.87 Abs growth (km 2 / yr) 0.61 0.02 0.01 0.04 0.03 Rel growth (%) 204 5 1 2 0 PHR=5 Area (km 2 ) 0.27 0.31 0.36 0.48 0.55 Abs growth (km 2 / yr) 0.27 0.02 0.05 0.05 0.05 Rel growth (%) 89 0 1 2 1 DEN=5 Area (km 2 ) 0.18 0.23 0.18 0.13 0.11 Abs growth (km 2 / yr) 0.18 0.02 0.04 0.03 0.01 Rel growth (%) 62 0 3 3 0 DEN=9 Area (km 2 ) 0.06 0.04 0.06 0.1 0.11 Abs growth (km 2 / yr) 0.06 0 0.02 0.02 0.01 Rel growth (%) 17 1 1 1 0 SURV=.02 Area (km 2 ) 0.46 0.5 0.6 0.8 0.92 Abs growth (km 2 / yr) 0.46 0.02 0.07 0.09 0.07 Rel growth (%) 150 1 2 2 1 SURV=.005 Area (km 2 ) 0.43 0.49 0.51 0.55 0.61 Abs growth (km 2 / yr) 0.43 0.02 0.01 0.02 0.03 Rel growth (%) 145 2 0 2 1 MAR=.65 Area (km 2 ) 0.04 0.08 0.08 0.08 0.08 Abs growth (km 2 / yr) 0.04 0.01 0 0 0.01 Rel growth (%) 15 0 0 1 0 MAR=.80 Area (km 2 ) 0.08 0.11 0.09 0.14 0.17 Abs growth (km 2 / yr) 0.08 0 0 0.03 0.01 Rel growth (%) 27 1 0 0 1
105 Table A 1. Continued SCR=High Profile Area (km 2 ) 0.24 0.26 0.82 1.33 1.86 Abs growth (km 2 / yr) 0.24 0.01 0.54 0.19 0.4 Rel growth (%) 79 1 25 4 9 DIST=100 Area (km 2 ) 0.44 0.66 0.73 0.82 0.93 Abs growth (km 2 / yr) 0.44 0.08 0.04 0.04 0.05 Rel growth (%) 148 3 1 1 1 DIST=300 Area (km 2 ) 0.43 0.6 0.7 0.85 0.99 Abs growth (km 2 / yr) 0.43 0.05 0.05 0.07 0.08 Rel growth (%) 143 0 0 1 1 Monte Carlo average over 100 simulation runs
106 APPENDIX B SUPPLEMENTARY DATA Figure B 1 E mpirical histograms for the number of new colonies formed by N reproductives (alates). Histograms are derived from 1000 model replications Gaussian probability distributions obtained using empirical means and standard d eviations are shown in yellow. In red the empirical kernel density estimates.
107 Figure B 1. Continued
108 Figure B 1. Continued
109 Figure B 1. Continued
110 Figure B 1. Continued
111 APPENDIX C COMPARISON OF MODEL OUTCOMES Table C 1 Simulation results over a homogeneous landscape. The average occupied areas and their standard deviations over 100 model replications are shown for the cell based and the IBM approaches, using different values for th e number of initial colonies N. Units are in hectares. Time N = 1 N = 5 N = 10 N = 25 N = 50 IBM CA IBM CA IBM CA IBM CA IBM CA Area SD Area SD Area SD Area SD Area SD Area SD Area SD Area SD Area SD Area SD 1 1 0 1 0 5 0 5.00 0 10 0 10.00 0 25 0 25.00 0 49 0 49.00 0 2 1.31 0.51 1 0 20.38 2.7 13.85 2.69 34.08 3.21 26.88 3.57 96.67 5.44 81.34 6.49 144.37 6.45 128.41 8.05 3 1.54 0.61 1.03 0.17 34.58 3.31 27.55 4.56 40.39 3.5 33.02 4.18 158.62 7.93 151.89 9.06 193.23 8.04 183.39 10.49 4 1.8 0.72 1.04 0.2 49.92 4.05 46.05 5.54 45.56 3.86 39.37 4.09 197.44 8.27 212.58 10.21 229.9 9.46 225.04 12.24 5 2.04 0.82 1.06 0.24 67.44 5.35 67.77 6.35 53.48 4.21 45.68 4.47 220.55 8.9 246.90 10.96 257.93 10.36 263.87 11.92 6 8.27 1.74 5.42 1.86 76.24 5.44 82.10 6.76 59.41 4.13 52.88 5.25 238.89 8.96 276.68 10.98 269.73 10.65 280.23 12.68 7 11.63 1.89 8.76 2.46 85.07 5.6 96.05 6.76 70.48 4.75 63.43 5.84 264.32 9.73 312.11 12.33 305.37 11.71 318.97 13.45 8 14 2.13 11.58 2.68 89.74 6.2 103.78 7.38 81.28 5.21 73.59 6.44 277.69 9.88 331.46 12.83 357.4 12.65 384.69 15.08 9 15.84 2.21 14.47 2.77 92.69 6.39 108.25 7.38 95.82 5.45 87.76 7.37 293.02 9.81 344.72 13.23 393.85 12.93 434.85 15.2 10 17.57 2.39 17.29 3.08 93.94 6.38 110.31 7.39 106.21 5.51 100.81 7.9 303.89 10.16 358.51 13.51 417.78 13.31 471.68 16.26 11 19.4 2.59 19.8 3.08 93.55 6.42 109.52 7.4 112.2 5.59 109.25 8.44 318.64 10.5 373.40 14.39 444.48 14.15 511.05 16.76 12 18.53 2.65 18.8 3.08 98.18 7.19 111.22 7.68 127.45 6.89 123.82 9.54 354.3 11.35 406.27 15.77 545.96 16.19 625.60 20.27 13 18.7 2.71 18.82 3.07 106.91 7.54 116.91 8.67 146.54 8.04 140.95 11.65 409.57 15.24 471.84 19.13 648.72 17.78 768.65 22.37 14 18.92 2.72 18.88 3.1 119.26 8.02 129.96 10.68 165.52 8.61 161.39 13.11 473.9 18.2 571.46 23.08 734.07 18.69 898.49 23.03 15 19.18 2.79 18.91 3.12 135.69 9.94 151.08 12.02 186.11 9.5 182.88 14.44 546.85 18.61 686.07 27.74 819.11 20.68 1020.14 23.94 16 20.93 3.07 19.49 3.39 152.06 11.13 177.51 14.69 204.02 10.9 206.39 16.11 608.08 20.09 793.93 29.39 904.73 20.94 1138.30 25.2 17 24.58 4.11 21.14 4.25 216.88 16.24 260.76 25.22 287.16 16.05 336.42 28.74 840.92 26.13 1101.08 36.59 1144.1 24.6 1444.79 36.97 18 29.16 5.25 23.78 5.39 278.18 17.61 362.51 30.07 344.05 17.04 434.50 34.38 1059.9 27.19 1402.46 35.41 1313.46 24.56 1629.30 37.47 19 33.88 6.23 27.08 6.65 341.96 18.37 476.39 33.01 392.85 19.15 514.23 34.82 1207.41 26.65 1590.05 32.6 1443.02 25.54 1766.70 37.36 20 38.28 7.01 32.05 7.98 403.34 18.85 586.05 34.14 453.34 20.75 612.53 35.55 1310.28 26.6 1717.78 31.46 1536.38 26.27 1872.97 34.82
112 APPENDIX D PRINCIPAL COMPONENTS LOADINGS Table D 1 Environmental and climate variables used in the study and principal component analysis (PCA) l oadings for all six extracted main axes, accounting for about 80% of the total variation in the original dataset. V ariables of the six axes are shown in bold. Variable Axis 1 Axis 2 Axis 3 Axis 4 Axis 5 Axis 6 X (Easting) 0.84 0.30 0.04 0.10 0.01 0.16 Y (Northing) 0.96 0.02 0.05 0.06 0.08 0.03 dew 0.97 0.11 0.06 0.01 0.06 0.08 prec 0.29 0.70 0.09 0.29 0.12 0.06 tmax 0.90 0.25 0.09 0.06 0.07 0.06 tmin 0.96 0.14 0.08 0.01 0.05 0.08 bio5 0.09 0.81 0.12 0.07 0.06 0.08 bio6 0.97 0.11 0.03 0.03 0.06 0.05 water 0.15 0.07 0.22 0.74 0.34 0.29 devel 0.19 0.33 0.75 0.11 0.16 0.33 barren 0.04 0.14 0.28 0.28 0.64 0.48 forest 0.76 0.04 0.06 0.01 0.15 0.13 shrub 0.42 0.30 0.00 0.06 0.27 0.49 herb 0.35 0.19 0.09 0.15 0.59 0.25 cultiv 0.17 0.56 0.13 0.58 0.02 0.48 wetlands 0.39 0.31 0.80 0.15 0.06 0.01 Eigenvalue 6.47 2.03 1.41 1.13 1.04 1.03 (%) Variance 40.43 12.66 8.81 7.06 6.48 6.42 Cumulative (%) Variance 40.43 53.09 61.90 68.96 75.44 81.86
1 13 APPENDIX E AST MODEL COEFFICIENTS Table E 1 N indicates the sample size used in each independent model replication. k indicates the number of i ndependent replications Only the best overall models for Florida are shown here. Model k N Parameter Coefficient SD LCI* UCI* BPOD MPC3 500 325 Intercept 16.67 1.12 18.79 14.5 1 Axis1 2.07 0.2 7 1.54 2.6 3 Axis2 0.81 0.44 0.02 1.7 2 Axis3 2.0 2 0.28 1.4 6 2.62 Prevalence 0.0 1 0.0 1 0 0.0 2 Maxent M1 16 10,000 Lambdas Axis1 14.48 Axis2 2.50 Axis3 10.16 *LCI = Lower Credible Interval (Bayesian lower confidence interval). UCI = Upper Credible Interval (Bayesian upper confidence interval)
114 APPENDIX F FST MODEL COEFFICIENTS Table F 1 F ST. N indicates the sample size used in each independent model replication. k indicates the number of i ndependent replications Only the best overall models for Florida are shown here. Mod el k N Parameter Coefficient SD LCI* UCI* BPOD MPC6 500 800 Intercept 5.60 0.52 6.8 7 4.83 Axis1 0.1 8 0.04 0.10 0.28 Axis2 0.75 0. 1 0.58 0.96 Axis3 1.24 0.2 4 0.93 1.90 Axis4 0.3 2 0.1 1 0.12 0.5 4 Axis5 0.4 5 0.1 1 0.27 0.67 Axis6 0.8 3 0. 20 0.55 1.3 5 Prevalence 0.0 5 0.02 0.01 0.10 Maxent MPC3 16 10,000 Lambdas Axis1 3.15 Axis2 3. 40 Axis3 6.1 3 *LCI = Lower Credible Interval (Bayesian lower confidence interval). UCI = Upper Credible Interval (Bayesian upper confidence interval)
115 LIST OF REFERENCES Adams, E.S., Levings, S.C., 1987. Territory size and population limits in mangrove termites. J. Anim. Ecol. 56, 1069 1081. Ajelli, M., Goncalves, B., Balcan, D., Colizza, V., Hu, H., Ramasco, J.J., et al., 2010. Comparing large scale computational approaches to epidemic modeling: agent based versus structured metapopulation models. BMC Infect. Dis. 10, 190. Akaike, H., 1974. A new look at the statistical model identification. IEEE Transactio ns on Automatic Control 19, 716 723. Alexanderian, A., Gobbert, M.K., Fister, K.R., Gaff, H., Lenhart, S., Schaefer, E., 2011. An age structured model for the spread of epidemic cholera: a nalysis and simulation. Nonlinear Anal. Real World Appl. 12, 3483 3498. Bordereau, C., Pasteels, J.M., 2011. Pheromones and chemical ecology of dispersal and foraging in termites. In: Bignell, D.E., Roisin, Y., Lo, N. (eds.), Biology of termites: a mod ern synthesis. Springer, pp. 279 320. Brown, D.G., Riolo, R., Robinson, D.T., North, M., Rand, W., 2005. Spatial process and data models: t oward integration of agent based models and GIS. J. Geogr. Syst. 7, 25 47. Burnham, B.O., 1973. Markov intertempor al land use simulation model. South. J. Agric. Econ. July, 253 258. Burnham, K.P., Anderson, D.R., 2002. Model s election and m ultimodel i nference: a p ractical i nformation t heoretic a pproach, 2nd ed. Springer Verlag. Carrasco, L.R., Mumford, J.D., MacLeo d, A., Harwood, T., Grabenweger, G., Leach, A.W., et al., 2010. Unveiling human assisted dispersal mechanisms in invasive alien insects: i ntegration of spatial stochastic simulation and phenology models. Ecol. Model. 221, 2068 2075. Chakraborty, A., Gelf and, A.E., Wilson, A.M., Latimer, A.M., Silander, J.A., 2011. Point pattern modelling for degraded presence only data over large regions. J. R. Stat. Soc.: Series C (App. Stat.) 5, 757 776. Clark, J.S., 1998. Why trees migrate so fast: confronting theory with dispersal biology and the paleorecord. Am. Nat. 152, 204 224. Clarke, K.C., Brass, J.A., Riggan, P.J., 1994. A cellular automaton model of wildfire propagation and extinction. Photogramm. Eng. Remote Sens. 60, 1355 1367.
116 Clarke, K.C., Gaydos, L.J ., 1998. Loose coupling a cellular automaton model and GIS: long term urban growth prediction for San Francisco and Washington/Baltimore. Int. J. Geogr. Inf. Sci. 12, 699 714. Cohen, J., 1960. A coefficient of agreement for nominal scales. Educ. Psychol. Meas. 20, 37 46. Collins, N.M., 1981. Populations, age, structure and survivorship of colonies of Macrotermes bellicosus (Isoptera: Macrotermitinae). J. Anim. Ecol. 50, 293 311. Conntable, S., Robert, A., Bordereau, C., 2012. Dispersal flight and colony development in the fungus growing termites Pseudacanthotermes spiniger and P. militaris Insect. Soc. 59, 269 277. Daly, C., Gibson, W.P., Taylor, G.H., Johnson, G.L., Pasteris, P., 2002. A knowledge based approach to the statistical mapping of climate. Clim. Res. 22, 99 113. Darlington, J.P.E.C., 1986. Seasonality in mature nests of the termite Macrotermes michaelse ni in Kenya. Insect. Soc. 33, 168 189. Della Pietra, S., Della Pietra, V., Lafferty, J., 1997. Inducing features of random fields. IEEE Trans. Pattern Anal. Mach. Intell. 19, 380 393. Diniz Filho, J.A., Bini, L.M., Rangel, T.F., Loyola, R.D., Hof, C., Nogus Bravo, D., et al., 2009. Partitioning and mapping uncertainties in ensembles of forecasts of species turnover under climate change. Ecogr. 32, 897 906. Divino, F., Golini, N., Jona Lasinio, G., Penttinen, A., 2013. Bayesian m odeling and MCMC c ompu tation in l inear l ogistic r egression for p resence only d ata. Cornell University Library. Source: http://arxiv.org/abs/1305.1232. Accessed 08/21/2013. Divino, F., Golini, N., Lasinio, G.J., Penttinen, A., 2011. Data a ugmentation a pproach in Bayesian m odell ing of p resence only d ata. Procedia Environ. Sci. 7, 38 43. Edwards, R., Mill, A.E., 1986. Termites in buildings: their biology and control. The Rentokil Library, 261 pp. Elith, J., Graham, C.H., Anderson, R.P., Dudk, M., Ferrier, S., Guisan, A., et al ., 2006. Novel methods improve prediction of species' distributions from occurrence data. Ecogr. 29, 129 151. Elith, J., Leathwick, J., 2007. Predicting species distributions from museum and herbarium records using multiresponse models fitted with multiv ariate adaptive regression splines. Divers. Distrib. 13, 265 275.
117 Elith, J., Leathwick, J.R., 2009. Species distribution models: ecological explanation and prediction across space and time. Annu. Rev. Ecol. Evol. Syst. 40, 677 697. Elith, J., Phillips, S.J., Hastie, T., Dudk, M., Chee, Y.E., Yates, C.J., 2011. A statistical explanation of MaxEnt for ecologists. Divers. Distrib. 17, 43 57. Evans, T.A., 2011. Invasive termites. In: Bignell, D.E., Roisin, Y., Lo, N. (eds.), Biology of termites: a modern synthesis. Springer, pp. 519 562. Farrar, D.E., Glauber, R.R., 1967. Multicollinearity in r egression a nalysis: t he p roblem r evisited. Rev. Econ. Stat. 49, 92 107. Fisher, R.A., 1937. The wave of advance of advantageous genes. Ann. Eugen. 7, 353 369. Franklin, J., 2009. Mapping species distributions: spatial inference and prediction. Cambridge University Press, Cambridge, UK. Gautam, B.K., Henderson, G., 2011. Relative h umidity p reference and s urvival of s tarved Formosan s ubterranean t ermites (Isopter a: Rhinotermitidae) at v arious t emperature and r elative h umidity c onditions. Env. Entomol. 40, 1232 1238. Grimm, V., Berger, U., Bastiansen, F., Eliassen, S., Ginot, V., Giske, J., et al., 2006. A standard protocol for describing individual based and agent based models. Ecol. Model. 198, 115 126. Grimm, V., Berger, U., DeAngelis, D.L., Polhill, J.G., Giske, J., Railsback, S.F., 2010. The ODD protocol: a review and first update. Ecol. Model. 221, 2760 2768. Grimm, V., Railsback, S.F., 2005. Individual based modeling and ecology. Princeton University Press, Princeton, NJ. Guisan, A., Thuiller, W., 2005. Pred icting species distribution: offering more than simple habitat models. Ecol. Lett. 8, 993 1009. Hanley, J.A., McNeil, B.J., 1982. The meaning and use of the area under a receiver operating characteristic (ROC) curve. Radiol. 143, 29 36. Hijmans, R.J., Cameron, S.E., Parra, J.L., Jones, P.G., Jarvis, A., 2005. Very high resolution interpolated climate surfaces for global land areas. Int. J. Clim. 25, 1965 1978. Hochmair, H.H., Scheffrahn, R.H., 2010. Spatial association of marine dockage with land born e infestations of invasive termites (Isoptera: Rhinotermitidae: Coptotermes ) in urban south Florida. J. Econ. Entomol. 103, 1338 1346.
118 Holmes, E.E., Lewis, M.A., Banks, J.E., Veit, R.R., 1994. Partial differential equations in ecology: spatial interaction s and population dynamics. Ecol. 75, 17 29. Hu, J., Zhong, J. H., Guo, M. F., 2007. Alate dispersal distances of the black winged subterranean termite Odontotermes formosanus (Isopera: Termitidae) in Southern China. Sociobiology 50 Husseneder, C., Simms, D.M., Ring, D.R., 2006. Genetic diversity and genotypic differentiation between the sexes in swarm aggregations decrease inbreeding in the Formosan subterranean termite. Insect. Soc. 53, 212 219. Huston, M., DeAngelis, D.L., Post, W.M., 1988. New computer models unify ecological theory. Biosci. 38, 682 691. Ikehara, S., 1966. Research report. Bull. Arts Sci. Div., 49 178. In Japanese. Jones, S.C., La Fage, J.P., Howard, R.W., 1988. Isopteran sex ratios: phylogenetic trends. Sociobiology 14, 89 1 56. Keller, L., 1998. Queen lifespan and colony characteristics in ants and termites. Insect. Soc. 45, 235 246. Lancaster, T., Imbens, G., 1996. Case control studies with contaminated controls. J. Econ. 71, 145 160. Law, A.M., Kelton, W.D., 1982. Sim ulation modelling and analysis. McGraw Hill Book Company, New York. Lee, S.H., Bardunias, P., Su, N.Y., 2008. Two strategies for optimizing the food encounter rate of termite tunnels simulated by a lattice model. Ecol. Model. 213, 381 388. Lelis, K., 20 06. Broward County surface water. Water Resources Division, Broward County Department of Natural Resource Protection. Source: http://gis.broward.org/GISData.htm. Accessed 05/25/2012. Leong, K.L.H., Tamashiro, Y.J., Su, N. Y., 1983. Microenvironmental fact ors regulating the flight of Coptotermes formosanus Shiraki in Hawaii (Isoptera: Rhinotermitidae). Proc. Hawaiian Entomol. Soc. 24 Lett, C., Silber, C., Barret, N., 1999. Comparison of a cellular automata network and an individual based model for the simu lation of forest dynamics. Ecol. Model. 121, 277 293.
119 Leuthold, R.H., Bruinsma, O., 1977. Pairing behavior in Hodotermes mossambicus (Isoptera). Psyche 84 Li, H. F., Fujisaki, I., Su, N. Y., 2013. Predicting h abitat s uitability of Coptotermes gestroi (I soptera: Rhinotermitidae) w ith s pecies d istribution m odels. J. Econ. Entomol. 106, 311 321. Li, H. F., Yang, R.L., Su, N.Y., 2010. Interspecific competition and territory defense mechanisms of Coptotermes formosanus and Coptotermes gestroi (Isoptera: Rhinotermitidae). Env. Entomol. 39, 1601 1607. Li, H. F., Ye, W., Su, N. Y., Kanzaki, N., 2009. Phylogeography of Coptotermes gestroi and Coptotermes formosanus (Isoptera: Rhinotermitidae) in Taiwan. Ann. Entomol. Soc. Am. 102, 684 693. Maidment, D.R., 1996. Environmental modeling within GIS. In: Goodchild, M.F., Steyaert, L.T., Parks, B.O., Johnston, C., Maidment, D., Crane, M., Glendinning, S. (eds.), GIS and environmental modelling: progress and research issues, Fort Collins, CO, pp. 3 15 324. Martius, C., 2003. Rainfall and air humidity: non linear relationships with termite swarming in Amazonia. Amazon. 17, 387 397. Messenger, M.T., Mullins, A.J., 2005. New flight distance recorded for Coptotermes formosanus (Isoptera: Rhinotermitid ae). Fla. Entomol. 88, 99 100. Mikler, A.R., Venkatachalam, S., Abbas, K., 2005. Modeling infectious diseases using global stochastic cellular automata. J. Biol. Syst. 13, 421 439. Mill, A.E., 1983. Observations on Brazilian termite alate swarms and so me structures used in the dispersal of reproductives (Isoptera: Termitidae). J. Nat. Hist. 17, 309 320. Morales Ramos, J.A., Rojas, M.G., 2005. Wood consumption rates of Coptotermes formosanus (Isoptera: Rhinotermitidae): a three year study using groups of workers and soldiers. Sociobiol. 45, 707 719. Morley, P.D., Chang, J., 2004. Critical behavior in cellular automata animal disease transmission model. Int. J. Mod. Phys. C 15, 149 162. Neubert, M.G., Kot, M., Lewis, M.A., 1995. Dispersal and pattern formation in a discrete time predator prey model. Theor. Popul. Biol. 48, 7 43. Noirot, C., 1990. C astes and reproductive strategies in termites. In: Engels, W. (ed.), Social insects: an evolutionary approach to castes and reproduction. Springer.
120 Nutti ng, W.L., 1969. Flight and colony foundation. In: Krishna, K., Weesner, F.M. (eds.), Biology of Termites, vol. 1. Academic Press, New York, pp. 233 282. Okubo, A., 1980. Diffusion and ecological problems: mathematical models. Biomath. 10 Peterson, A.T., Ortega Huerta, M.A., Bartley, J., Sanchez Cordero, V., Soberon, J., Buddemeier, R.H., et al., 2002. Future projections for Mexican faunas under global climate change scenarios. Nat. 416, 626 629. Peterson, C.J., 2010. Termites: here, there and everywhere? EARTH magazine 55, 46. Phillips, S.J., Anderson, R.P., Schapire, R.E., 2006. Maximum entropy modeling of species geographic distributions. Ecol. Model. 190, 231 259. Phillips, S.J., Dudik, M., 2008. Modeling of species distributions with Maxent: new extensions and a comprehensive evaluation. Ecogr. 31, 161 175. Phillips, S.J., Dudk, M., Schapire, R.E., 2004. A maximum entropy approach to species distribution modeling. Proc. 21st Int. Conf. M ach. Learn. 655 662 Pitt, J.P.W., 2008. Modelling the spread of invasive species across heterogeneous landscapes. Lincoln University, p. 232. Pitt, J.P.W., 2009. Predicting Argentine ant spread over the heterogeneous landscape using a spatially explicit stochastic model. Ecol. Appl. 19, 1176 1186. Pitt, J.P.W., Kriticos, D.J., Dodd, M.B., 2011. Temporal limits to simulating the future spread pattern of invasive species: Buddleja davidii in Europe and New Zealand. Ecol. Model. 222, 1880 1887. R Develop ment Core Team, 2013. A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria. Ramirez Villegas, J., Jarvis, A., 2010. Downscaling g lobal c irculation m odel o utputs: t he d elta m ethod. Decis. Policy Anal Working Pap. No. 1 Robert, C.P., Casella, G., 2004. Monte Carlo s tatistical m ethods, 2nd ed. Springer Verlag New York, Inc., Secaucus, NJ. Rust, M.K., Su, N. Y., 2012. Managing s ocial i nsects of u rban i mportance. Annu. Rev. Entomol. 57, 355 375. Scheffrahn, R.H., 2013. Overview and c urrent s tatus of n on n ative t ermites (Isoptera) in Florida. Fla. Entomol. 96, 781 788.
121 Scheffrahn, R.H., Busey, P., Edwards, J.K., Krecek, J., Maharajh, B., Su, N. Y., 2001. Chemical prevention of colony foundation by Cryptotermes brevis (Isoptera: Kalotermitidae) in attic modules. J. Econ. Entomol. 94, 915 919. Scheffrahn, R.H., Cabrera, B.J., Kern Jr, W.H., Su, N. Y., 2002. Nasutitermes costalis (Isoptera: Termitidae) in Florida: first record of a non endemic estab lishment by a higher termite. Fla. Entomol. 85, 273 275. Scheffrahn, R.H., Crowe, W., 2011. Ship borne termite (Isoptera) border interceptions in Australia and onboard infestations in Florida, 1986 2009. Fla. Entomol. 94, 57 63. Scheffrahn, R.H., Hochmair, H.H., Kern Jr, W.H., Warner, J., Krecek, J., Maharajh, B., et al. 2014. T argeted elimination of the exotic termite, Nasutitermes corniger (Isoptera: Termitidae: Nasutitermitinae), from infested tracts in southeastern Florida Int J. Pest Manag. Source: http://dx.doi.org/10.1080/09670874.2014.882528 Accessed 02/28/2014 Scheffrahn, R.H., Su, N.Y., 2005. Distribution of the termite genus C optotermes (Isoptera: Rhitotermitidae) in Florida. Fla. Entomol. 88, 201 203. Skellam, J.G., 1951. Random dispersal in theoretical populations. Biom. 38, 196 218. Steyaert, L.T., 1993. A perspective for studying of environmental simulation. In: Goodchild, M.F., Parks, B.O., Steyaert, L.T. (eds.), Environmental modelling with GIS. Oxford University Press, New York, pp. 16 30. Su, N. Y., Scheffrahn, R.H., 1987. Alate production of a field colony of the Formosan subterranean termite (Isoptera: Rhinotermitidae). Sociobiology 13, 209 215. Su, N. Y., Scheffrahn, R.H., 2000. Formosan subterranean termite, Coptote rmes formosanus Shiraki (Insecta: Isoptera: Rhinotermitidae). Institute of Food and Agricultural Sciences, University of Florida. Source: http://edis.ifas.ufl.edu/in278. Accessed 05/23/2012. Su, N.Y., Scheffrahn, R.H., 1997. A new introduction of a subter ranean termite, Coptotermes havilandi Holmgren (Isoptera: Rhinotermitidae) in Miami, Florida. Fla. Entomol. 80, 408 411. Sudshira, H.S., 2010. Integration of agent based and cellular automata models for simulating urban sprawl. International Institute fo r Geo Information Science and Earth Observation, p. 78. Syfert, M.M., Smith, M.J., Coomes, D.A., 2013. The effect of sampling bias and model complexity on the predictive performance of maxent species distribution models. PLoS ONE 8
122 Tanner, M., Wong, W., 1987. The calculation of posterior distribution by data augmentation. J. Am. Stat. Assoc. 82, 528 550. Thorne, B.L., 1982. Termite termite interactions: workers as an agonistic caste. Psyche 89, 133 150. Thorne, B.L., 1983. Alate production and sex rat io in colonies of the neotropical termite Nasutiternes corniger (Isoptera; Termitidae). Oecol. 58, 103 109. Toffoli, T., Margolus, N., 1987. Cellular automata machines: a new environment for modelling. MIT Press, Cambridge, MA. Tonini, F., Hochmair, H.H., Scheffrahn, R.H., DeAngelis, D., 2013. Simulating the spread of an invasive termite in an urban environment using a stochastic individual based model. Environ. Entomol. 42, 412 423. University of Florida GeoPlan Center, 2010. Generalized land use derived from 2010 parcels Florida DOT District 4. Source: http://www.fgdl.org. Accessed 01/24/2014 van den Bosch, E., Metz, J.A.J., Diekmann, O., 1990. The velocity of spatial population expansion. J. Math. Biol. 28, 529 565. Ward G., Hastie, T., Barry, S.C., Elith, J., Leathwick, J.R., 2009. Presence only data and the EM algorithm. Biom. 65, 554 563. Warren, D.L., Seifert, S., 2011. Environmental niche modeling in Maxent: the importance of model complexity and the performance o f model selection criteria. Ecol. Appl. 21, 335 342. only data in ecology. Ann. Appl. Stat. 4, 1383 1402. Wiegand, T., Revilla, E., Knauer, F., 2004. Dealing with uncertainty in spatially explicit population models. Biodivers. Conserv. 13, 53 78.
123 BIOGRAPHICAL SKETCH Francesco Tonini was born in Florence, Italy. In 2007, h e graduated with a Bachelor o f Science in s tatistics and i nformation s ystems from the University of Florence, Italy. In 2010, he graduated with a Master of Science in s tatistics for biomedicine, the environment, and technology from the University of Rome "La Sapienza", Italy. During t he same year, he was awarded a four year Graduate School Fellowship Program by the University of Florida and accepted as a Ph.D. student. The award funded Francesco's position as a graduate research assistant between the years 2010 2014. In 2014, he graduated with a Ph.D. in f orest r esources and c onservation with a concentration in g eomatics
xml version 1.0 encoding UTF-8
REPORT xmlns http:www.fcla.edudlsmddaitss xmlns:xsi http:www.w3.org2001XMLSchema-instance xsi:schemaLocation http:www.fcla.edudlsmddaitssdaitssReport.xsd
INGEST IEID EZNGG7CTU_L2KJKA INGEST_TIME 2014-10-03T21:29:40Z PACKAGE UFE0046461_00001
AGREEMENT_INFO ACCOUNT UF PROJECT UFDC
Ecological Modelling 276 (2014) 1Â… 13 Contents lists available at ScienceDirect Ecological Modelling jo ur nal ho me page: www.elsevier.com/locate/ecolmodel Impacts of deer management practices on the spatial dynamics of the tick Ixodes ricinus : A scenario analysisSen Li a , Sophie O. Vanwambeke a Alain M. Licoppe b Niko Speybroeck c aGeorges Lematre Centre for Earth and Climate Research, Earth and Life Institute, Universit catholique de Louvain, Louvain-la-Neuve, BelgiumbDepartment of Natural and Agricultural Environment Studies, Operational Directorate-General for Agriculture, Natural Resources and the Environment (DGARNE), Gembloux, BelgiumcInstitute of Health and Society (IRSS), Universit catholique de Louvain, Brussels, Belgium a r t i c l e i n f o Article history: Received 13 June 2013 Received in revised form 12 December 2013 Accepted 21 December 2013 Keywords: Cellular automata Deer management Spatial heterogeneity Tick population ecology Tick controla b s t r a c t Deer, for example roe deer, red deer and fallow deer, are the common reproduction host types for European Ixodes ricinus ticks. Understanding the consequences of deer management on the spatial dynamics of ticks may advise the risk management of tick-borne diseases, and thus be of public health importance. We present a scenario analysis to understand such consequences by integrating multi-disciplinary knowledge within a predictive modelling framework. A spatial tick population model was adopted to explore how changes in the host population may affect woodland patchand landscape-level tick dynamics. Scenarios on current and foreseen European deer management strategies were built based on expert knowledge. These scenarios were then tested through the described model for their potential effectiveness as tick control strategies. Our models indicate that: (i) reducing local deer densities could not effectively reduce tick abundance if woodland patches are well-connected, allowing deer population exchanges, (ii) controlling deer grazing intensity in grassland may not be an effective tick control strategy, (iii) local extinction of deer could decrease tick abundance considerably but deer reintroduction could lead to fast tick upsurge, and (iv) controlling human disturbances may reduce the tick density at landscape-level, as well as tick ÂhotspotsÂŽ (i.e., areas with unusually high tick density) at woodland patch-level. Our results can instruct policy-makers on the potential impact on public health of wildlife management strategies, and suggest empirical investigations of disease risks. For optimising such simulation studies on disease risks, however, a better understanding of how biodiversity may inÂ”uence the ecology of tick and pathogen transmission is required. 2014 Elsevier B.V. All rights reserved. 1. Introduction Ixodes ricinus (Acari: Ixodidae) is the most abundant and widely distributed tick species in Europe. Over the past two decades, the reported incidence of tick-borne diseases, such as Lyme borreliosis and tick-borne encephalitis, has increased in Europe ( Randolph, 2008 ). This has made tick-borne zoonoses the most prevalent vector-borne diseases in Europe and hence kindled public health and scientiÂ“c attention ( Morens et al., 2004 ). Host availability is fundamental for the maintenance of ticks and pathogen populations ( Ruiz-Fons and Gilbert, 2010 ). Success Corresponding author at: 457, Place Louis Pasteur 3, B1348 Louvain-la-Neuve, Belgium. Tel.: +32 0 10472870; fax: +32 0 10472877. E-mail addresses: email@example.com firstname.lastname@example.org (S. Li), email@example.com (S.O. Vanwambeke), firstname.lastname@example.org (A.M. Licoppe), email@example.com (N. Speybroeck). in managing tick-borne diseases depends on understanding tick dynamics, including the way in which tick feeding success is affected by varying host population densities and distributions. Ixodes ticks feed on a broad range of hosts. I. ricinus larvae primarily feed on small-sized animals such as rodents and insectivores ( Gray et al., 1994; Jaenson and Talleklint, 1996 ), while hosts for nymphs are less speciÂ“c, including birds, and small, medium, and large-sized mammals ( Estrada-Pe na et al., 2005 ). Adult ticks often have a narrower host range, with mediumand large-sized hosts, for example deer ( Rizzoli et al., 2009; Tagliapietra et al., 2011 ). This three-host tick life cycle can be impacted by a diapause occurring at different life stages such that its duration may extend over two years ( Gray, 1981 ). Ixodes ticks can be transported across the landscape as hosts move. Landscape effects, e.g., woodland fragmentation, shaping host distribution and limiting host movement at different spatial scales, can yield different inÂ”uences on tick populations in different life stages ( Li et al., 2012a ). Deer and other large mammals are the target of direct management practices (e.g., fencing, hunting, translocation, etc.). In previous empirical studies, 0304-3800/$ Â… see front matter 2014 Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.ecolmodel.2013.12.023
2 S. Li et al. / Ecological Modelling 276 (2014) 1 Â… 13 deer density has been found to be positively related to tick densities and the tick infestations on rodents ( Cagnacci et al., 2012; Killilea et al., 2008 ). There is a need to analyse the potential consequences of these management practices on the spatial dynamics of ticks, so that both wildlife managers and the public would be better informed about potential changes in disease risks ( Luckhart et al., 2010 ). Predicting the impact of deer management on the spatial dynamics of I. ricinus requires a thorough understanding of (i) deer management practices, (ii) the spatially explicit response of deer movement patterns to such management practices and (iii) the spatial ecology of tick-deer relations. Agent-based modelling and cellular automata have the capacity to integrate these aspects ( Grimm et al., 2005 ). To date, only a few studies have reported on biological process-based models modelling tick populations in a spatially explicit way ( Killilea et al., 2008 ). These models have been used mainly to explore how spatial dynamics of tick densities and/or infections can be inÂ”uenced by landscape heterogeneity and host movements. Some recent examples are: a multi-habitat model ( Hoch et al., 2010 ) for I. ricinus tick, an agent-based model ( Wang et al., 2012 ) for Amblyomma americanum a multi-patch model ( Watts et al., 2009 ) and a reaction diffusion model ( Jones et al., 2011 ) for louping-ill virus, and a cellular automata model for Lyme disease ( Li et al., 2012a ). In general, these studies support either a single host following a random exploratory movement or assume a certain number of hosts continuously travelling between woodland patches. However, in the Â“eld, animals can exhibit a set of statistically distinguishable movement patterns as a function of their internal state, motion capacity, navigation capacity, environmental conditions and the management strategies impacting them ( Nathan et al., 2008 ). Therefore, a multiphasic hypothesis for the host movement pattern (i.e., considering more than one host movement phase, e.g., both home-ranging phase, in which animals randomly move within their home range, and displacement phase, in which animals tend to move over longer distances for new habitats) may be an interesting extension to the current spatial tick modelling context. This study aimed to investigate the consequences of deer management practices on the spatial dynamics of ticks. We Â“rstly described a spatial tick model adapted from a previous study ( Li et al., 2012a ). Then, based on expert knowledge, four scenarios of current and foreseen deer management practices were built, namely: (i) reducing local deer density, (ii) controlling deer grazing intensity in grassland, (iii) translocation of deer species, and (iv) controlling human disturbance and deer displacement between woodland patches. The inÂ”uence on the spatial dynamics of ticks under these scenarios was tested through the model. Finally, we discussed these results with a focus on the effectiveness of these practices as tick controlling strategies. 2. Materials and methods 2.1. Model proÂ“le The cellular automata model for the spatial ecology of ticks was adapted from Li et al. (2012a) and included signiÂ“cant modiÂ“cations compared to the original model. Firstly, the tick feeding pattern was more detailed to better represent the reality by enabling low probabilities for adults to feed on small mammals ( Cagnacci et al., 2012 ) and larvae to feed on deer ( Kiffner et al., 2010 ). The simpliÂ“ed method based on Â“xed tick attachment rates in the original model was replaced by a more detailed one based on host-Â“nding probabilities of ticks and controlled by the feeding capacities of ticks on hosts. This adaption further allows to simulate the effects of deer density, as a tick ampliÂ“er, on tick infestation levels ( Cagnacci et al., 2012 ). Secondly, rules for deer movement were developed to distinguish behaviours in home-ranging and displacement phases, allowing to represent more realistic deer movement patterns. Finally, pathogen transmission functions were not included as the study at hand focused on tick populations. The modiÂ“cations aimed at improving the modelÂs speciÂ“city, as the focus in this study was on the effects of the movements of the major host types on the abundance and distribution of ticks. The model was coded into Repast Simphony ( North et al., 2006 ). The space used was two-dimensional and organised as a lattice of cells. Each cell in the lattice had three layers: a tick population, a host population and a landscape layer. The lattice was programmed to evolve in discrete time steps following a set of transition rules to update the cell states simultaneously. The present model adopts a cell size of 1 ha and a time step of 1 week. 2.1.1. States and parameters The three layers of cell states were: Tick population layer : A ÂlarvaÂŽÂ…ÂnymphÂŽÂ…ÂadultÂŽ life stage structure was used ( Fig. 1 ). In each stage, ticks could be in questing, feeding or interstadial development phases. When encountering hosts, questing ticks attached for blood meals, then dropped and developed into the next life stages. Female adult ticks (assuming half of the emerged adult ticks were females) also needed blood meals to produce (i.e., lay eggs that hatch into) larval ticks. Host population layer: Two generalised host types were used: small mammals (including rodents and lagomorphs) and deer. Each life stage of the ticks was assumed capable of feeding on both host types. The overall number of the two host types was Â“xed but the host distributions could vary between time steps as a result of movements. Indeed, host movements resulted in ticks being transported from one place to another. Landscape layer: The cell states used were ÂwoodlandÂŽÂ… ÂgrasslandÂŽÂ…Ânon-vegetated areasÂŽ. Deer were considered to inhabit woodland mostly and to spend a proportion of time in grassland. Ticks could inhabit grassland but woodland was considered more suitable. Small mammals could inhabit both woodland and grassland. Mortality rates of ticks, densities of hosts and movement patterns of hosts varied with land cover types. Hosts could enter non-vegetated areas, but would not stay, meaning that no ticks would drop off in non-vegetated area. The landscape layer was set static and could be based on either real or theoretical landscape maps. Parameters in the model were estimated from Â“eld and laboratory observations found in the literature ( Table 1 ). 2.1.2. Transition rules 184.108.40.206. Rules for modelling the tick population development. For each cell at time step t we modelled the questing tick populations as (see Table 1 for the values used and the sources of parameters): qLt= (1 ÂŠ mqL) qLt ÂŠ 1+ 0 5 (1 ÂŠ mAL)dAL (1 ÂŠ mfL) ( fAD t ÂŠ dALÂŠ 1+ fAS t ÂŠ dALÂŠ 1) ÂŠ fLD t ÂŠ 1ÂŠ fLS t ÂŠ 1(1) qNt= (1 ÂŠ mqN) qNt ÂŠ 1+ (1 ÂŠ mLN)dLN (1 ÂŠ mfN) ( fLD t ÂŠ dLNÂŠ 1+ fLS t ÂŠ dLNÂŠ 1) ÂŠ fND t ÂŠ 1ÂŠ fNS t ÂŠ 1(2) qAt= (1 ÂŠ mqA) qAt ÂŠ 1+ (1 ÂŠ mNA)dNA (1 ÂŠ mfA) ( fND t ÂŠ dNAÂŠ 1+ fNS t ÂŠ dNAÂŠ 1) ÂŠ fAD t ÂŠ 1ÂŠ fAS t ÂŠ 1(3)
S. Li et al. / Ecological Modelling 276 (2014) 1 Â… 13 3 Fig. 1. Tick life cycle in the model ModiÂ“ed from Li et al. (2012a) where qL qN and qA indicate the questing populations of larval, nymphal and adult ticks; mqL, mqNand mqAthe weekly tick mortality survival rates in questing phases; mfL, mfNand mfAthe weekly tick mortality survival rates in feeding phases; dAL, dLNand dNAare development periods between feeding and moulting into the next life stages and mAL, mLNand mNAare weekly tick mortality rates in those periods. All m and d parameters are constant. All mortality rates are considered four times higher in grassland than in woodland ( Sf ), in relation to the drier conditions in grassland. indicates the number of eggs per adult. The factor of 0.5 in Eq. (1) indicates that half of the adult ticks are assumed to be eggs-producing females. fL fN and fA indicate whether the feeding populations were larval, nymphal and adult ticks. Superscripts D and S indicate the tick population is feeding on deer and small mammals. The number of weekly feeding ticks in each cell was estimated by fTH= qTt pTH, where qT refers to the questing tick population of a life stage ( L N A ) and pTH, the host-Â“nding probabilities, refers to the probability of ticks in the life stage ( T ) to attach on a host type ( D S ). The estimated feeding ticks per host were assumed not to exceed the feeding capacity on the host ( C ), or the maximum number of tick attachments on one host in one week. Such capacities were estimated based on Â“eld observations of the ranges of tick attachments per host and tick feeding durations ( Table 1 ). For example, given that the duration of the larva attachment can be as short as two days and the number of larvae on one small mammal can be as many as 30 ( Cagnacci et al., 2012; Gray, 2002 ), the weekly feeding capacity of small mammals for larvae ( CLS) was estimated as (30 larvae/2 days) 7 days = 105 larvae. Feeding durations were assumed to be less than one week for larvae (approximately three days, c.f. Macleod (1932) ), and one week for nymphs and adults ( Hancock et al., 2011 ). 220.127.116.11. Rules for modelling the host movement and tick transportation in a cellular space. A multi-phasic movement pattern was considered: (i) in home-ranging phase, animals randomly move within their home range, while (ii) in displacement phase, animals tend to move over longer distances for new habitats. Such a multiphasic hypothesis has been utilised in both theoretical models ( Bartumeus et al., 2005; Skalski and Gilliam, 2003; Zhang et al., 2007 ) and empirical studies on movement ecology ( Fryxell et al., 2008; Morales et al., 2004 ). Deer, for example roe deer Capreolus capreolus red deer Cervus elaphus and fallow deer Dama dama are capable of moving over long distances. In the model, deer populations were initially considered in home-ranging phase but a proportion were allowed to switch to a displacement phase ( pDis ) when competition or human disturbance increased ( Nathan et al., 2008 ). During the home-ranging phase, deer were allowed to explore within and around woodland. They could also venture into grassland for short stays ( Putman, 1986 ), i.e., a proportion of a time step ( pG ). It has been reported that, the home ranges of roe deer, red deer and fallow deer can be as large as 1 km2, 7.5 km2or 20 km2( Putman, 1988 ). During the displacement phase, deer were allowed to move over long distances (i.e., up to 5 km for roe deer and 120 km for red deer per week has been reported) looking for new woodland habitats that are not over-populated (or below the limit of the carrying capacity, K of the woodland cell) ( Georgii and Schroder, 1983; Mysterud, 1999; Wahlstrom and Liberg, 1995 ). In this phase, deer were assumed to travel quickly through grassland and the time spent in grassland was considered negligible. An extended Moore neighbourhood ( Fig. 2 ) was considered for host movement. The host movement capacity ( MC ) denotes the minimum distance travel from the centre to the boundary of the neighbourhood. Thus the neighbourhood sizes were different for different host types in different movement phases. Two sets of rules were designed, respectively, for host movement patterns in home ranging and displacement phases. Movement rules in home ranging phase : The extended neighbourhood in this phase represents the home range of the host type concerned. By assuming that hosts had an equal probability to
4 S. Li et al. / Ecological Modelling 276 (2014) 1 Â… 13 Table 1 Model parameters. Parameter Symbol Where a Value Source Parameters in relation to tick life cycle Average no. eggs per adult W/G 2000 Randolph and Craine (1995) Mortality rate in woodland (weekÂŠ 1) Questing larvae mqLW 0.03 Daniel et al. (1976) Questing nymphs mqNW 0.03 Daniel et al. (1976) Questing adult ticks mqAW 0.02 Daniel et al. (1976) Developing from engorged larvae into questing nymphs mLNW 0.03 Daniel et al. (1976) Gray (1981) Developing from engorged nymphs into questing adults mNAW 0.01 Daniel et al. (1976) Gray (1981) Developing from engorged adults into questing larvae mALW 0.02 Daniel et al. (1976) Gray (1981) Scaling factor for mortality rates of questing and developing ticks Sf G 4 Mount et al. (1997) Talleklint and Jaenson (1997) Feeding mortality rate (weekÂŠ 1) Larvae mfLH 0.5 Gray (1981) Hancock et al. (2011) Nymphs mfNH 0.5 Gray (1981) Hancock et al. (2011) Adults mfAH 0.35 Gray (1981) Hancock et al. (2011) Duration of developmental phases (week) Developing from engorged larvae into questing nymphs dLNW/G 45 Macleod (1932) Developing from engorged nymphs into questing adults dNAW/G 53 Macleod (1932) Developing from engorged adults into questing larvae dALW/G 45 Macleod (1932) Small mammals Â“nding probability Questing larvae pLSW/G 0.3 Hancock et al. (2011) Questing nymphs pNSW/G 0.2 Hancock et al. (2011) Questing adults pASW/G 0.01 Assumption Deer Â“nding probability Questing larvae pLDW/G 0.001 Assumption Questing nymphs pNDW/G 0.1 Hancock et al. (2011) Questing adults pADW/G 0.3 Hancock et al. (2011) Host feeding capacity (weekÂŠ 1) Maximum larva attachments on one small mammal CLSH 105 Cagnacci et al. (2012) Gray (2002) Maximum nymph attachments on one small mammal CNSH 6 Oliver (1989) Cagnacci et al. (2012) Maximum adult attachments on one small mammal CASH 0.1 Cagnacci et al. (2012) Oliver (1989) Maximum larva attachments on one deer CLDH 800 Gray (2002) Kiffner et al. (2010) Maximum nymph attachments on one deer CNDH 300 Oliver (1989) Kiffner et al. (2010) Maximum adult attachments on one deer CADH 150 Oliver (1989) Kiffner et al. (2010) Parameters in relation to host movement patterns Movement capacity (m weekÂŠ 1) Small mammal in home ranging phase MCSW/G 100 Kikkawa (1964) Deer in home ranging phase MCD HW 500 Putman (1988) Deer in displacement phase MCD DW 2000 Georgii and Schroder (1983) Mysterud (1999) Wahlstrom and Liberg (1995) Proportion of time step spent in grassland for deer (%) pG G 35 Putman (1986) Proportion of deer population in displacement phase (%) pDis W 0.1 Assumption Parameters in relation to landscape heterogeneity Density of small mammal (haÂŠ 1) dSW/G 80 Escutenaire et al. (2000) Density of deer (haÂŠ 1) dDW 0.12 Delbeuck (2008) OFFH Carrying capacity of small mammal (haÂŠ 1) KSW/G 120 Assumption Carrying capacity of deer (haÂŠ 1) KDW 0.3 Assumption a This column indicates where the parameter can be applied: W, in woodland; G, in grassland, H, on host animals. move in each direction, the home ranging population of a host type in the original cell was evenly divided into eight parts to be considered for movement. In each time step, the following two sub-steps (i) and (ii), respectively, for the trigger and completion of movement were repeated eight times, once for each direction. In sub-step (i), a location of the destination was randomly generated by using a modiÂ“ed Gaussian model (a detailed description can be found in Appendix S2 in Li et al. (2012a) ) giving the concerned home ranging host population (i.e., the one eighth part) a higher probability to move at shorter distances. A movement trigger value of 100 m (size of one sell) was adopted, representing the required distance to reach the adjacent cells. Thus, the hosts would stay in
S. Li et al. / Ecological Modelling 276 (2014) 1 Â… 13 5 Fig. 2. Extended Moore neighbourhood. There are eight directions: east, northeast, north, northwest, west, southwest, south, and southeast. the original cell if the distance to the randomly-generated destination was less than the movement trigger value. If the distance was higher than the hostsÂ movement capacity (with a low probability), the destination would be adjusted to the closest boundary cell of the extended neighbourhood, so that the movements beyond home ranges were disallowed. In sub-step (ii), there were two different situations when considering the completion of movement. First, completion of movement was possible when the destination was a suitable habitat for the concerned host type. The hosts could then only complete the movement if the host population was below the carrying capacity of its host type in the destination cell. Second, completion of movement was not possible when the destination was a less suitable habitat for the concerned host type (i.e., when considering the deer moving to grassland). The hosts would then spend a proportion of time in the cell and then return to the original cell. In all other situations, the movements in this direction were considered to have failed, with as a result hosts remaining in the original cell. Movement rules in displacement phase : Host populations in displacement phase were considered to move together from a cell. The extended neighbourhood in this phase represents the search window of the host. Displacing populations would start to search from the nearest eight neighbouring cells to the next 16 adjacent cells and then reach over the extended neighbourhood. The searching activity would be Â“nished when a habitat cell was found in a different patch and the host population there was below the host carrying capacity of the cell. If no such habitat cell was found, the hosts remained in the original cell at the end of the time step. As we focused on deer management, small mammals were assumed to perform home ranging behaviours only. It has been reported that small mammals have relatively small home range areas (e.g., about 0.2 ha for the bank vole Myodes glareolus and the wood mice Apodemus sylvaticus ( Kikkawa, 1964; Kozakiewicz et al., 2007 ). In the model, small mammals were assumed to inhabit both woodland and grassland and their movement capacity ( MCS) was assumed to be 100 m. Given a value of 100 m to trigger the movement between adjacent cells, it was therefore modelled that most small mammal movements could be completed within each cell, except for a small proportion, amongst those inhabiting areas close to cell boundaries. Moving hosts transported ticks. As the time step of this study was one week, nymphs and adults, assumed to take one week to feed, could be transported between cells. The transport of larvae, however, was considered negligible, as larvae were assumed to feed for less than one week. When a proportion of hosts had moved between habitats, the same proportion of feeding ticks was transported. When deer had only spent a proportion of a time step in grassland, the same proportion of ticks feeding on them dropped off, and the same proportion of questing ticks were picked up. 2.1.3. Process overview and scheduling At each successive time step, transition rules for tick development and host movement patterns were applied. For each life stage, questing ticks that emerged (moulted in the previous time step from a previous life stage) were added in. Questing ticks that died in the previous time step were removed. Then, ticks attached on hosts in the same cell and host movements were considered. Attached ticks on the out-moving animals were transported and all dropped off at the end of the time step. Cell states were updated simultaneously at each time step. 2.1.4. Inputs, outputs and initialisation The model required the input of a land cover map and outputted the questing population of ticks in different life stages. The land cover map could be a real world map or an artiÂ“cial map designed for theoretical investigations. In the present study, all simulations were initialised with a tick density of 50 000 adults per ha in woodland areas, according to a range of 0Â…1620 adults sampled per ha in Belgium ( Li et al., 2012b ) and by assuming a 5Â…9% sampling efÂ“ciency ( Daniels et al., 2000 ). Hosts were assumed homogenously distributed in habitats at the initial time step. Then, as the model ran and movement rules were applied, their populations would become heterogeneous and varying between time steps. 2.2. Model evaluation 2.2.1. Comparison with Â“eld observations The model performance was examined by testing its ability to reÂ”ect relative differences in observed tick abundance in different sites in south Belgium as reported in Li et al. (2012b) Tick sampling data for comparison were selected (in total 125 samples in 51 sites) based on the following criteria: (i) in Wallonia (deer data); (ii) in rural areas (to avoid potential strong effects of human disturbance on local deer density); and (iii) sampled in the same period of the same year (to avoid potential strong inÂ”uence of climatic differences on tick population). Using these criteria, only four samples could be used (Thuin, Marche-en-Famenne, Wellin and Marcheles-Dames). Land cover maps (5000 m 5000 m, resolution 100 m) included ÂwoodlandÂŽ (including broadleaf and coniferous forests), ÂgrasslandÂŽ (including pasture and moorland) and Ânon-vegetated areaÂŽ (including built-up areas, agricultural land and forest clearings). Deer densities in the woodland were calculated based on the spring standing population of roe, red and fallow deer provided by the ÂDirection Gnrale des Ressources Naturelles et de lÂEnvironnementÂŽ of the Walloon Region. 2.2.2. Sensitivity analysis The sensitivity of the model to all parameters was assessed with the land cover map of Thuin (i.e., the site with the greatest number of sampled ticks) as used in Section 2.2.1 The sensitivity index was calculated as: S = Log10( Di/ D0)/Log10( Pi/ P0), with D0the density index at equilibrium when using the default values for all parameters (values in Table 1 ) and Dithe density index when the value of a parameter ( P ) is increased from the default P0to Piby 5%. It has been demonstrated in a number of biological models ( Keeling and Gilligan, 2000; Ogden et al., 2005 ) that S is a robust indicator of the sensitivity. A higher absolute value of S means a stronger effect of the parameter change on the tick abundance. S values of 1 and ÂŠ 1 indicate positive and negative linear effects.
6 S. Li et al. / Ecological Modelling 276 (2014) 1 Â… 13 Fig. 3. ConÂ“guration of the artiÂ“cial landscape. The spatial extension of the landscape is 50 50 cells (cell size = 1 ha). Grey area indicates three woodland patches: (A) size = 47 ha, (B) size = 177 ha, and (C) size = 226. The minimum movement capacities required for deer to move between A and B is 800 m, between B and C is 500 m, and between A and C is 500 m. White areas indicate grassland. 2.3. Scenarios Four scenarios relevant to current and foreseen deer management practices in Europe were developed. SpeciÂ“cally, the scenarios address the potential effects of wildlife management practices on deer patterns. The model predicted their consequences on the spatial pattern of ticks using an artiÂ“cial landscape ( Fig. 3 ). It further allows assessing whether I. ricinus can be managed through controlling its major reproduction hosts. The design of the theoretical landscape was based on the research questions to be explored. As the main target of the present study was the spatial relationship between tick and deer populations, the structure of the theoretical landscape was kept simple to limit the effects of landscape conÂ“guration on the results. The woodland area was assumed to be fragmented, reÂ”ecting a common Â“eld situation. Grassland and woodland covers were included, as these are the major vegetated land cover types in many European countries. Such design allowed testing the effects of deer movement between-land cover types and between woodland patches. 2.3.1. Scenario 1 Â… reducing local deer density Reducing the deer density, via hunting, can sometimes beneÂ“t conservational plans for the protection of woodland resources and for a better recovering and quality of deer population. In the simulation, we intended to examine whether a failed deer density reduction in one woodland patch will inÂ”uence the tick density over the landscape. Deer densities observed by the ÂDirection Gnrale des Ressources Naturelles et de lÂEnvironnementÂŽ in Belgium were used to parameterise the model. We assumed the three woodland patches had a deer density ( dD) of 0.18 heads/ha (i.e., close to the highest deer density of 0.19 heads/ha in Bllingen in 2008) and the controlling target was 0.08 heads/ha (i.e., the averaged deer density in South Belgium in 2008). Then, four situations were considered: (i) the control was successful in all woodland patches, and a situation for which the controlling failed (i.e., deer density did not change) (ii) in woodland A, (iii) in woodland B or (iv) in woodland C. For each situation, three sets of simulations were compared with movement capacity of deer ( MCD H) = 800 m (i.e., deer can move between woodlands A, B and C), 500 m (i.e., deer can move between woodlands A and C and between woodlands B and C but not between woodlands A and B) and 300 m (i.e., deer cannot reach to other woodland patches). Our hypothesis was that, if deer can move between woodland patches to transport ticks, the failed reduction of deer in one woodland patch may increase the density in both local and nearby woodland patches, resulting in a general increase in tick density at landscape level. 2.3.2. Scenario 2 Â… controlling deer grazing intensity in grassland Grassland grazing by deer, especially by red and fallow deer, has been reported to cause signiÂ“cant damage on agricultural grassland production in Europe ( Licoppe, 2006; Putman and Moore, 1998; Trdan and Vidrih, 2008 ). Roe deer, the most abundant and selective browser of European ungulate species, is often underestimated in its impact on grassland, pasture and ecotone composites of trees and shrubs ( Putman and Moore, 1998; Reimoser and Putman, 2011 ). The prevention of such damage can be reached via exclusive fencing. As complete removal of deer grazing in grassland by fencing is considered undesirable in many natural forests, there is an essential need to understand the impact of grazing intensities on woodland dynamics for potential deer management policies ( Hester et al., 2000 ). In reality, managing the grazing intensity of wildlife is not easy, but via hunting the numbers can be regulated, allowing for density changes in speciÂ“c areas ( Hester et al., 2000 ). Moreover, a cline of grassland usage may also reÂ”ect the situations in which enclosures are incomplete. Thus, a proportion of deer are free-ranging deer that have not been enclosed, or the deer fenced inside the woodland can jump over the low fence and still graze in grassland. In the model, grazing intensity was represented in two ways: the grazing range and time spent in grassland. In the simulation, the impact of reduced deer grazing intensities on tick dynamics was tested by decreasing movement capacities of deer ( MCD H= 400 m, 300 m, 200 m and 100 m) and the proportion of time spent in grassland ( pG = 75%, 50%, 25% and 0% of a week). Thus, in all simulations, deer populations could not exchange between woodland patches. Our hypothesis was that controlled deer grazing intensity could lead deer to spend more time in woodland and thus increase the tick density in woodland. 2.3.3. Scenario 3 Â… translocation of deer species Translocation of deer species has been widely practised in Europe during the recent past for meeting the local hunting demands and rehabilitating habitats where local deer species were extinct. Roe deer have been reintroduced in some of their historic regions in France ( Calenge et al., 2005 ), Italy ( Dupanloup et al., 2002 ) and Portugal ( Vernesi et al., 2002 ). The current distribution of a western red deer lineage in eastern Europe was suggested to be mainly artiÂ“cial and a result of translocation ( NiedziaÂ›kowska et al., 2010 ). In the model, we tested the effects of different time periods between the local extinction and reintroduction of deer species. Deer populations in all woodlands were removed after 10 years of simulation. Then reintroductions were considered after 1, 2, 3, 5 and 10 years, respectively. The aim was to explore to what extent the deer extirpation could reduce the local tick population and how long it would take to restore the local tick population after deer reintroduction. Our hypothesis was that the extinction of deer species would not eradicate ticks (as tick population could still be maintained by small mammals), but would greatly reduce the tick population (as only limited numbers of adult ticks could feed on small mammals). 2.3.4. Scenario 4 Â… controlling human disturbance and deer displacement between woodland patches Human disturbance (e.g. hunting, tourism, management practice etc.) may result in deer displacement and a net deer population exchange between woodland patches. In the displacement phase, red deer may migrate to other woodland patches, hide for several days and then return ( Sunde et al., 2009 ). Roe deer, though preferring to hide in ground vegetation, may also perform displacement when being released after capture for Â“tting tracking devices ( Morellet et al., 2009 ). Besides, immigration of other deer species, for example via natal dispersal ( Prvot and Licoppe, 2013 ), may result in the competitive displacement of roe deer ( Dolman and Wber, 2008; Latham, 1999 ). In the model, controlling human disturbance was assumed to reduce deer displacement rates. We assumed up to 50% of deer population to be displaced ( pDis ) by disturbance ( Sunde et al., 2009 ) and a lower movement capacity ( MCD H= 300 m) for home-ranging deer to avoid disturbance
S. Li et al. / Ecological Modelling 276 (2014) 1 Â… 13 7 Table 2 Validation results. Site Deer density (haÂŠ 1) No. of sampled nymphs Estimated nymph density in Â“eld (haÂŠ 1) a Simulated nymph density (haÂŠ 1) Thuin 0.13 125 2.78Â…5.00 E4 2.95 0.01 E4 Marche-en-Famenne 0.14 108 2.40Â…4.32 E4 3.00 0.01 E4 Wellin 0.11 88 1.96Â…3.52 E4 2.41 0.01 E4 Marche-les-Dames 0.07 42 0.93Â…1.68 E4 1.52 0.00 E4 a Assuming 5Â…9% sampling efÂ“ciency ( Daniels et al., 2000 ). ( Jayakody et al., 2011 ). Finally, we simulated a set of decreased deer displacement rates ( pDis from 50% to 0% of deer population in displacement phase, with a successive decrease of 5%) for testing their impact on tick population dynamics. Our hypothesis was that a reduced displacement rate of deer would reduce the tick density, as it could relatively increase the home-ranging population of deer (so that more ticks would be transported to grassland in which tick mortalities are higher). 3. Results Equilibrium levels (with a weekly change of questing tick densities <0.5%) were achieved in all simulations in Â“ve to eight simulated years after initialisation. Hence, to have a safe conÂ“rmation of the equilibrium, outcomes were recorded from the 10th to 30th simulated years (521stÂ…1560th weeks) for model evaluation, sensitivity analysis and scenario analysis. Different initial values of tick densities in different life stages were tested as well and did not affect the results. 3.1. Model evaluation 3.1.1. Comparison with Â“eld observations The results in Table 2 show that, the simulated output largely reÂ”ects the differences of Â“eld sampled data in the four sites (comparing Â“eld sampled data and simulated data: PearsonÂs chisquared = 6.05, p -value = 0.1). The unexplained difference (e.g., the sample from Thuin shows the highest numbers of ticks while in Marche-en-Famenne numbers are highest with the model) may be explained by different reasons. Microclimatic conditions, for example, may have been different on the sampling occasions and local small mammal densities can vary. The local deer compositions may also have inÂ”uenced the results, as, based on the Â“eld data, Thuin and Marche-les-Dames were frequented by roe deer only, while Marche-en-Famenne and Wellin were frequented by both roe and red deer. 3.1.2. Sensitivity analysis In general, tick densities in different post-egg life stages were sensitive to different sets of parameters ( Fig. 4 ). Larval tick densities were sensitive to average numbers of eggs per adult, the weekly mortality rate of questing larvae, the duration of and the weekly mortality rate in the developmental phase from engorged adult ticks to questing larvae, the density of deer, and the weekly maximum number of adult attachments on one deer. Nymphal tick densities were sensitive to many parameters. High sensitivities were observed for the mortality in developing phases from engorged larvae to questing nymphs and from engorged adults to questing larvae, the mortality rate of feeding larvae, the mortality rate of questing nymphs, the weekly maximum number of larvae and nymphs attachments on one small mammal, and the density of deer and small mammals. Adult tick densities show high sensitivities to developing phase mortality from engorged nymphs to questing adults. They were also sensitive to the weekly mortality rates in questing adult phase, the weekly maximum nymph attachment on small mammals and small mammal densities. 3.2. Scenario analysis 3.2.1. Scenario 1 Â… reduction of deer population may not be an effective method to control the tick population if deer can move between neighbouring woodland patches A successful deer population reduction from 0.18 to 0.08 heads/ha in all woodland patches contributed to a reduction of approximately 50.1% of the tick population at landscape level ( Fig. 5 a and e). When deer were able to move between woodland patches (movement capacity = 800 m and 500 m), failing to reduce deer density in one woodland patch would increase the questing tick density in all other woodland patches ( Fig. 5 bÂ…d). There was no inÂ”uence on tick densities in other woodlands, when deer could not move to neighbouring patches (movement capacity = 300 m). Overall questing tick densities were highest when the reduction of deer population failed in woodland C. Woodland C, being relatively larger, would hold a larger deer population for the same deer density. No major changes in tick abundance in grassland were observed. 3.2.2. Scenario 2 Â… controlled deer grazing intensity can reduce the tick abundance in grassland, but increase it in woodland As deer grazing intensity was controlled (or deer movement capacity and the proportion of time steps deer spent in grassland decreased) in the simulations, questing tick abundance decreased in grassland but increased in woodland ( Fig. 6 ). Extensive deer browsing control in grassland (by decreasing the movement capacity from 400 to 100 m per week and the proportion of time steps spent in grassland from 75% to 0%) could reduce the questing tick abundance in grassland by 87.1% and increase it by 10.7% in woodland, resulting in a combined overall landscape tick increase effect of 10.5% (the simulated ratio of tick density in grassland to that in woodland was approximately 1:105). 3.2.3. Scenario 3 Â… translocation of deer: local extinction of deer could decrease tick abundance, but could not eradicate ticks Local deer extinction could decrease the questing tick density substantially, namely by 76% (failing to Â“nd hosts before dying) ( Fig. 7 ). There were two equilibrium levels: (i) a high level of 3.5 105questing ticks per ha with both deer and small mammal populations and (ii) a low level of 8.5 104questing ticks per ha when only small mammals were present. In all situations, the local questing tick abundance started to decrease after 45 weeks of deer extirpation and to increase 45 weeks after deer reintroduction. This is associated with the model parameters on the duration of developmental phases ( Table 1 ). When tick started to decrease after deer extirpation, it took about 140 weeks to reduce to a lower equilibrium level. When tick started to increase after deer reintroduction, it took about 210 weeks to recover the tick density to a higher equilibrium level. 3.2.4. Scenario 4 Â… decreased displacement of deer could decrease the tick density in woodland, and also avoid the formation of areas with extremely high tick density Decreased deer displacement decreased average questing tick densities in woodland patches, as well as the highest tick density
8 S. Li et al. / Ecological Modelling 276 (2014) 1 Â… 13 Fig. 4. Sensitivity of model variables. Sensitivity indices were calculated for woodland larval tick density ( Sl ), nymphal tick density ( Sn ) and adult tick density ( Sa ). in woodland cells ( Fig. 8 ). A control of displacement proportion of deer population from 50% to 0% reduced the overall questing tick density from 3.9 105to 3.6 105ticks per ha and reduced the highest questing tick density from 7.4 106to 4.6 105ticks per ha. Along with an increased displacement rate, fewer ticks were noticed in grassland (not shown). 4. Discussion This study aimed to explore possible consequences of wildlife management practices on tick spatial dynamics. A populationbased, multi-host and multi-land cover cellular automata model for the I. ricinus population ecology was adapted from a previous
S. Li et al. / Ecological Modelling 276 (2014) 1 Â… 13 9 Fig. 5. Simulated questing tick densities under different deer population reduction scenarios. MC = movement capacity in home ranging phase. Failed control means the deer population remained at 0.18 heads/ha. The control target was 0.08 heads/ha. study ( Li et al., 2012a ) so that it allowed simulating how various deer movement patterns inÂ”uence the spatial distribution and abundance of ticks in heterogeneous landscapes. A set of scenarios on current and foreseen deer management practices and the corresponding changes of deer movement patterns in Europe were established. Those scenarios were then tested for their potential impacts on the spatial dynamics of ticks through the model using an artiÂ“cial landscape, comprised of three woodland patches surrounded by grassland. We showed how complex ecological models for disease vectors can be adapted for exploring consequences on tick densities of various management strategies, allowing policy makers to assess the potential effects of managing deer population on public health. The study Â“rstly predicted that controlling the deer population would not always be effective to control ticks. In isolated areas, reducing the deer population can reduce the local tick abundance considerably even if the deer reduction attempts failed in other areas. However, when target areas were well-connected with other areas for deer population exchange, the reduction of local tick abundance would not be effective unless deer populations were reduced in all connected areas. In line with a number of Â“eld studies ( Estrada-Pe na, 2002, 2003 ), this prediction emphasises the effect of interactions between landscape connectivity and spatial heterogeneity of host populations. Moreover, the relative location of woodland patches in the network appeared to be inÂ”uential on the spatial tick population patterns. In scenario 1, while deer could not directly move between woodland A and B (movement capacity = 500 m), the failed control of deer population in woodland A still increased the tick population in woodland B via woodland C. This may suggest a direction for future empirical work towards examining the potential of geographical accessibility of woodland patches (the degree to which the woodland is available to as many deer as possible) as a local risk factor for ticks and tick-borne diseases. To help preventing tick-borne diseases, hunting, for example, may be planned as a means to control deer population in woodland of great geographical accessibility. A second prediction was that decreasing deer grazing in grassland adjacent to forest moderately reduces the tick abundance in grassland and increases the tick abundance in woodland. This is in accordance with previous modelling studies ( Hoch et al., 2010; Jones et al., 2011 ), suggesting a higher tick density in enclosed forests ( Gilbert et al., 2012 ). Ticks increased at the landscape level, as the exposure of ticks to the higher mortality rates in grassland decreased. Moreover, ticks feeding on rodents may be increased where deer presence is strictly controlled (e.g., via fencing), which may result in an increase in disease transmission ( Perkins et al., 2006; Ros and Pugliese, 2007 ). Hence, controlling deer browsing between connecting and contrasting land covers may not be an effective way to reduce tick abundance and the risk of exposure to pathogens. Besides, along with previous Â“ndings from scenario 1, the two game-harvesting strategies in Europe could be compared: hunting free ranging animals or animals in enclosed forest parks. Hunting in parks is prohibited in Belgium but is allowed in some countries such as France, Spain and Portugal. The access to pastures or other woodlands is impossible for park game animal and the tick abundance could be higher in parks than in free ranging
10 S. Li et al. / Ecological Modelling 276 (2014) 1 Â… 13 Fig. 6. Relative abundance of questing ticks in grassland (a) and in woodland (b) under different deer browsing controlling scenarios. Deer grazing scenarios were set up by using all possible combinations of deer movement capacity (400 m, 300 m, 200 m and 100 m) and proportion of time deer spent in grassland (75%, 50%, 25% and 0%). The relative tick abundance was calculated by dividing the tick abundance simulated for each scenario by the default tick abundance simulated (i.e., when both movement capacity and proportion of time spent in grassland were the highest). areas. Promoting empirical tick investigations in forest parks may be necessary for better acknowledging potential contacts between ticks and humans. A third prediction was that the local extinction of deer species could lead to an approximately 76% control of tick populations in Fig. 8. Effects of deer displacement on the questing tick densities in woodland. Averaged density ( ) and maximum density ( ) were shown. woodland. This is slightly lower than the empirical 86Â…96% Â“eld control of ticks in Scottish forests and moorland by excluding deer locally, through fencing ( Gilbert et al., 2012 ). The observed difference may be due to stochasticity but also factors such as a higher tick mortality rates ( Macleod, 1932 ) or a relative lower small mammal abundance ( Gilbert et al., 2000 ) found in moorland may have played a role here. In our simulations, the tick population was maintained by an unchanged number of smaller mammals across the landscape. However, there is evidence that the exclusion of deer species in woodland can result in an increase of small mammals while the introduction of deer species may not affect local small mammal communities ( Smit et al., 2001 ). Moreover, the sensitivity analysis indicated that an increased small mammal population can increase the nymphal and adult tick abundance ( Fig. 4 ). Reintroducing deer may raise tick abundance to a higher level than before. Hence, the reintroduction of deer species into tick-infested areas should be designed carefully, especially concerning the location of the release sites. A fourth prediction was that a decreased deer displacement between forest patches in response to a control of human disturbance could decrease the possibility of tick ÂhotspotsÂŽ with particular high tick densities, as well as the overall tick densities. Controlling human disturbance of animals thus seemed to be effective in controlling ticks. In scenario 4, deer exchange between woodland patches resulted from deer displacement only, as homeranging deer could not reach other woodland patches. Increased deer displacement could reduce the proportion of home-ranging Fig. 7. Evolutions of questing tick densities in woodland under different deer translocation scenarios. Deer were extirpated in the simulated week of 521 (10 years). 1 year, 2 years, 3 years, 5 years and 10 years indicate deer reintroductions after 1, 2, 3, 5 and 10 years of deer extirpation.
S. Li et al. / Ecological Modelling 276 (2014) 1 Â… 13 11 deer, leading to fewer ticks being transported from woodland to adjacent grassland where tick mortalities were high. Hence, the effect of deer displacement on tick abundance was positive at landscape-level. The positive local-level effect on the maximum tick density was due to the simulated clustered distribution of deer in woodland, which is in accordance with Â“eld observations under high human disturbance ( Jayakody et al., 2008 ). In the model, all deer followed the same rules to select destination woodland cells for displacement. Some cells were selected constantly by moving deer as they are closer (and thus easier) to travel to. This resulted in deer population clusters and therefore more drop-offs of fed ticks in and around these cells. As a result, ÂhotspotsÂŽ of ticks were formed. In general, deer displacement seemed to not play an important role at landscape level, which is in line with the Â“ndings in the sensitivity analysis. However, its effects on the local tick distribution could be vast and deems further exploration. Such Â“ndings could hardly be gained from traditional non-spatial models, highlighting the utility of spatial ecological models in the present study. Finally, there is a public health reason to create Âquiet zonesÂŽ (i.e., area excluding human disturbance) in a relatively isolated location to reduce local deer movement. Such management practice can lower the contact probability between animals and humans. In any case, the spatial heterogeneity found in and between results underlines the need to know those places effectively used by humans, where the latter will be exposed to tick bites. The model could be extended in different ways, allowing a better understanding e.g., of the impact of local biodiversity. Firstly, the host layer could be adapted in a species-speciÂ“c fashion. In the present study, a generalised layer mixing several common deer species was used to cope with current data limitations with respect to the composition of local deer populations. However, roe, red and fallow deer are different in size and in their capacities for hosting ticks ( Matuschka et al., 1993; Vor et al., 2010 ), and can establish different movement patterns across multiple spatio-temporal scales. In Wallonia, the authorities aim to reduce the red deer and wild boar densities (see legal rule M.B.08.05.1993 amended in 2008 by Decree of the Walloon Government), which could have an indirect positive effect on roe deer populations. It may be interesting to study how such modiÂ“cation of deer composition could inÂ”uence tick spatial dynamics. Secondly, the risk of tick-borne diseases, such as Lyme disease and tick-borne encephalitis, could be modelled by including a disease transmission function allowing exploring the tick infection prevalence. I. ricinus ticks can feed on a wide range of host species. However, some of these hosts are incompetent for pathogen transmission, such as deer, whilst most small mammals, especially rodents and squirrels, are competent vectors. Therefore, deer are commonly hypothesised to have a Âdilution effectÂŽ, thus an increase in deer population can decrease the probability of tick feeding on small mammals and hereby decrease the tick infection prevalence. However, such a hypothesis has been reported to be more complex than explained and scale-dependent ( Cagnacci et al., 2012 ). Future spatial modelling studies may allow exploring whether or not the ÂdilutionÂŽ effects depend on landscape conÂ“gurations. Improved understanding of the spatio-temporal dynamics of tick feeding patterns on hosts is also desirable. Still, such developments of the model rely heavily on interdisciplinary investigations, for example between ecology, geography and computer sciences, on the local composition and movement patterns of wildlife species. Thirdly, future modelling attempts may focus on Â“ner spatial and temporal scales. A more detailed tick biology model could be daily-based, so that the feeding durations and the tick transportation during feeding can be represented in a better way. On a Â“ner spatial scale, more landscape effects can be explored, for example the presence and characteristics of ecotones (i.e., woodland-grassland interfaces) that are often favourable to tick survival and development. Finally, it can be promising to include the phenology of tick. The seasonal population dynamics of ticks and their host animals have been studied extensively ( Killilea et al., 2008; Tagliapietra et al., 2011 ). A number of tick behaviours have been reported to be highly sensitive to season, for example diapause, maturation, and host-Â“nding activity. The used modelling approach has some shortcomings resulting from the simplifying assumptions made, for example about the tick biology and climatic inÂ”uence. This has also been reÂ”ected in the sensitivity analysis: the most sensitive variables concerned tick mortality in different phases (i.e., questing, feeding and developing between life stages) and feeding tick numbers on hosts. Similar model sensitivities have been found for existing models, e.g., the I. scapularis model of Ogden et al. (2005) and I. ricinus model of Hoch et al. (2010) Both models were highly sensitive to the feeding mortality of larvae and feeding numbers of immature ticks. For a better understanding of the spatial tick dynamics, a better representation of tick biology in natural conditions and an improved function for the development of tick life stages may be needed. To our knowledge, there is no clear evidence to claim that density dependence acts on I. ricinus numbers. The use of this factor in models can perfectly balance the populations but this can be equally well done through climate-related mortalities only. Andrewartha and Birch (1954) dismiss the existence of density dependence in the population dynamics of I. ricinus In our models we used density dependence, but future modelling may need to include climate to investigate whether or not climate is sufÂ“cient to explain the Â“eld observations. Further empirical investigations could help to increase the scientiÂ“c knowledge about the aforementioned aspects. However, a model is always a purposeful representation of reality. Including additional unknowns may increase the model complexity by creating too many degrees of freedom and thus can be harmful. Hence, there is a need in the design of ecological models to maximise the speciÂ“city while minimising the complexity. 5. Conclusions In summary, we used a spatial modelling approach to integrate existing interdisciplinary knowledge for predicting consequences of different wildlife management strategies. Prior information and expert knowledge on deer movement patterns and their possible responses to the management strategies were combined to improve the understanding of the spatial dynamics of deer, the tick natural hosts. Then, by using a cellular automata model, the effects of several potential wildlife management practices were evaluated such as the control of the deer population, grazing, translocation and displacement on the spatial dynamics of tick at woodland patchand landscape-levels. The modelling approach presented here can be used for further experiments about the adaptive management of wildlife and tick-borne diseases. Acknowledgements Sen Li is a research fellow (Aspirant FNRS) at the National Fund for ScientiÂ“c Research, Belgium. The authors wish to thank the anonymous reviewers for their valuable comments that have enhanced the work. ReferencesAndrewartha, H.G., Birch, L.C., 1954. The Distribution and Abundance of Animals. University of Chicago Press. Bartumeus, F., Da Luz, M.G.E., Viswanathan, G.M., Catalan, J., 2005. Animal search strategies: a quantitative random-walk analysis. Ecology 86, 3078Â…3087. Cagnacci, F., Bolzoni, L., Ros, R., Carpi, G., Hauffe, H.C., Valent, M., Tagliapietra, V., Kazimirova, M., Koci, J., Stanko, M., Lukan, M., Henttonen, H., Rizzoli, A., 2012. Effects of deer density on tick infestation of rodents and the hazard of tick-borne encephalitis, I: empirical assessment. Int. J. Parasitol. 42, 365Â…372.
12 S. Li et al. / Ecological Modelling 276 (2014) 1 Â… 13 Calenge, C., Maillard, D., Invernia, N., Gaudin, J.C., 2005. Reintroduction of roe deer Capreolus capreolus into a Mediterranean habitat: female mortality and dispersion. Wildl. Biol. 11, 153Â…161. Daniel, M., Cerny, V., Dusbabek, F., Honzakova, E., Olejnicek, J., 1976. InÂ”uence of microclimate on the life cycle of the common tick Ixodes ricinus (L.) in thermophilic oak forest. Folia Parasitol. 23, 327Â…342. Daniels, T.J., Falco, R.C., Fish, D., 2000. Estimating population size and drag sampling efÂ“ciency for the blacklegged tick (Acari: Ixodidae). J. Med. Entomol. 37, 357Â…363. Delbeuck, C., 2008. Cellule etat de lÂenvironnement Wallon: Environmental Outlook for Wallonia 2008. Service Public de Wallonie. Dolman, P.M., Wber, K., 2008. Ecosystem and competition impacts of introduced deer. Wildl. Res. 35, 202Â…214. Dupanloup, I., Schneider, S., ExcofÂ“er, L., 2002. A simulated annealing approach to deÂ“ne the genetic structure of populations. Mol. Ecol. 11, 2571Â…2581. Escutenaire, S., Chalon, P., Verhagen, R., Heyman, P., Thomas, I., Karelle-Bui, L., AvsicZupanc, T., Lundkvist, A., Plyusnin, A., Pastoret, P.P., 2000. Spatial and temporal dynamics of Puumala hantavirus infection in red bank vole ( Clethrionomys glareolus ) populations in Belgium. Virus Res., 91Â…107. Estrada-Pe na, A., 2002. Understanding the relationships between landscape connectivity and abundance of Ixodes ricinus Ticks. Exp. Appl. Acarol. 28, 239Â…248. Estrada-Pe na, A., 2003. The relationships between habitat topology, critical scales of connectivity and tick abundance Ixodes ricinus in a heterogeneous landscape in northern Spain. Ecography 26, 661Â…671. Estrada-Pe na, A., Osacar, J.J., Pichon, B., Gray, J.S., 2005. Hosts and pathogen detection for immature stages of Ixodes ricinus (Acari: Ixodidae) in North-Central Spain. Exp. Appl. Acarol. 37, 257Â…268. Fryxell, J.M., Hazell, M., Borger, L., Dalziel, B.D., Haydon, D.T., Morales, J.M., McIntosh, T., Rosatte, R.C., 2008. Multiple movement modes by large herbivores at multiple spatiotemporal scales. Proc. Natl. Acad. Sci. U.S.A. 105, 19114Â…19119. Georgii, B., Schroder, W., 1983. Home range and activity patterns of male red deer ( Cervus elaphus L.) in the Alps. Oecologia 58, 238Â…248. Gilbert, L., Jones, L.D., Hudson, P.J., Gould, E.A., Reid, H.W., 2000. Role of small mammals in the persistence of Louping-ill virus: Â“eld survey and tick co-feeding studies. Med. Vet. Entomol. 14, 277Â…282. Gilbert, L., Maffey, G.L., Ramsay, S.L., Hester, A.J., 2012. The effect of deer management on the abundance of Ixodes ricinus in Scotland. Ecol. Appl. 22, 658Â…667. Gray, J.S., 1981. The fecundity of Ixodes ricinus (L.) (Acarina, Ixodidae) and the mortality of its developmental stages under Â“eld conditions. Bull. Entomol. Res. 71, 533Â…542. Gray, J.S., 2002. Biology of Ixodes species ticks in relation to tick-borne zoonoses. Wien. Klin. Wochenschr. 114, 473Â…478. Gray, J.S., Kahl, O., Janetzkimittman, C., Stein, J., Guy, E., 1994. Acquisition of Borrelia burgdorferi by Ixodes ricinus ticks fed on the European hedgehog, Erinaceus europaeus L. Exp. Appl. Acarol. 18, 485Â…491. Grimm, V., Revilla, E., Berger, U., Jeltsch, F., Mooij, W.M., Railsback, S.F., Thulke, H.H., Weiner, J., Wiegand, T., DeAngelis, D.L., 2005. Pattern-oriented modeling of agent-based complex systems: lessons from ecology. Science 310, 987Â…991. Hancock, P.A., Brackley, R., Palmer, S.C., 2011. Modelling the effect of temperature variation on the seasonal dynamics of Ixodes ricinus tick populations. Int. J. Parasitol. 41, 513Â…522. Hester, A.J., Edenius, L., Buttenschn, R.M., Kuiters, A.T., 2000. Interactions between forests and herbivores: the role of controlled grazing experiments. Forestry 73, 381Â…391. Hoch, T., Monnet, Y., Agoulon, A., 2010. InÂ”uence of host migration between woodland and pasture on the population dynamics of the tick Ixodes ricinus : a modelling approach. Ecol. Model. 221, 1798Â…1806. Jaenson, T.G.T., Talleklint, L., 1996. Lyme borreliosis spirochetes in Ixodes ricinus (Acari: Ixodidae) and the varying hare on isolated islands in the Baltic Sea. J. Med. Entomol. 33, 339Â…343. Jayakody, S., Sibbald, A.M., Gordon, I.J., Lambin, X., 2008. Red deer Cervus elephus vigilance behaviour differs with habitat and type of human disturbance. Wildl. Biol. 14, 81Â…91. Jayakody, S., Sibbald, A.M., Mayes, R.W., Hooper, R.J., Gordon, I.J., Lambin, X., 2011. Effects of human disturbance on the diet composition of wild red deer ( Cervus elaphus ). Eur. J. Wildl. Res. 57, 939Â…948. Jones, E., Webb, S., Ruiz-Fons, F., Albon, S., Gilbert, L., 2011. The effect of landscape heterogeneity and host movement on a tick-borne pathogen. Theor. Ecol. 4, 435Â…448. Keeling, M.J., Gilligan, C.A., 2000. Bubonic plague: a metapopulation model of a zoonosis. Proc. R. Soc. Lond. B: Biol. 267, 2219Â…2230. Kiffner, C., Ldige, C., Alings, M., Vor, T., Rhe, F., 2010. Abundance estimation of Ixodes ticks (Acari: Ixodidae) on roe deer ( Capreolus capreolus ). Exp. Appl. Acarol. 52, 73Â…84. Kikkawa, J., 1964. Movement, activity and distribution of the small rodents Clethrionomys glareolus and Apodemus sylvaticus in woodland. J. Anim. Ecol. 33, 259Â…299. Killilea, M., Swei, A., Lane, R., Briggs, C., Ostfeld, R., 2008. Spatial dynamics of Lyme disease: a review. Ecohealth 5, 167Â…195. Kozakiewicz, M., Choluj, A., Kozakiewicz, A., 2007. Long-distance movements of individuals in a free-living bank vole population: an important element of male breeding strategy. Acta Theriol. 52, 339Â…348. Latham, J., 1999. InterspeciÂ“c interactions of ungulates in European forests: an overview. For. Ecol. Manag. 120, 13Â…21. Li, S., Hartemink, N., Speybroeck, N., Vanwambeke, S.O., 2012a. Consequences of landscape fragmentation on Lyme disease risk: a cellular automata approach. PLoS ONE 7, e39612. Li, S., Heyman, P., Cochez, C., Simons, L., Vanwambeke, S., 2012b. A multi-level analysis of the relationship between environmental factors and questing Ixodes ricinus dynamics in Belgium. Parasit. Vectors 5, 149. Licoppe, A., 2006. The diurnal habitat used by red deer ( Cervus elaphus L.) in the Haute Ardenne. Eur. J. Wildl. Res. 52, 164Â…170. Luckhart, S., Lindsay, S.W., James, A.A., Scott, T.W., 2010. Reframing critical needs in vector biology and management of vector-borne disease. PLoS Negl. Trop. Dis. 4, e566. Macleod, J., 1932. The bionomics of Ixodes ricinus L. the sheep tick of Scotland. Parasitology 24, 382Â…400. Matuschka, F.R., Heiler, M., Eiffert, H., Fischer, P., Lotter, H., Spielman, A., 1993. Diversionary role of hoofed game in the transmission of Lyme disease spirochetes. Am. J. Trop. Med. Hyg. 48, 693Â…699. Morales, J.M., Haydon, D.T., Frair, J., Holsiner, K.E., Fryxell, J.M., 2004. Extracting more out of relocation data: building movement models as mixtures of random walks. Ecology 85, 2436Â…2445. Morellet, N., Verheyden, H., Angibault, J.-M., Cargnelutti, B., Lourtet, B., Hewison, M.A.J., 2009. The effect of capture on ranging behaviour and activity of the European roe deer Capreolus capreolus Wildl. Biol. 15, 278Â…287. Morens, D.M., Folkers, G.K., Fauci, A.S., 2004. The challenge of emerging and reemerging infectious diseases. Nature 430, 242Â…249. Mount, G.A., Haile, D.G., Daniels, E., 1997. Simulation of blacklegged tick (Acari: Ixodidae) population dynamics and transmission of Borrelia burgdorferi J. Med. Entomol., 461Â…484. Mysterud, A., 1999. Seasonal migration pattern and home range of roe deer ( Capreolus capreolus ) in an altitudinal gradient in southern Norway. J. Zool. 247, 479Â…486. Nathan, R., Getz, W.M., Revilla, E., Holyoak, M., Kadmon, R., Saltz, D., Smouse, P.E., 2008. A movement ecology paradigm for unifying organismal movement research. Proc. Natl. Acad. Sci. U.S.A. 105, 19052Â…19059. NiedziaÂ›kowska, M., Jedrzejewska, B., Honnen, A.C., Otto, T., Sidorovich, V.E., Perzanowski, K., Skog, A., Hartl, G.B., Borowik, T., Bunevich, A.N., Lang, J., Zachos, F.E., 2010. Molecular biogeography of red deer Cervus elaphus from eastern Europe: insights from mitochondrial DNA sequences. Acta Theriol., 1Â…12. North, M.J., Collier, N.T., Vos, J.R., 2006. Experiences creating three implementations of the repast agent modeling toolkit. ACM Trans. Model. Comput. Simul. 16, 1Â…25. OFFH, LÂObservatoire de la faune, de la Â”ore et des habitats en Wallonie Ministre de la Rgion Wallonne, DGRNE. Ogden, N.H., Bigras-Poulin, M., OÂCallaghan, C.J., Barker, I.K., Lindsay, L.R., Maarouf, A., Smoyer-Tomic, K.E., Waltner-Toews, D., Charron, D., 2005. A dynamic population model to investigate effects of climate on geographic range and seasonality of the tick Ixodes scapularis Int. J. Parasitol., 375Â…389. Oliver Jr., J.H., 1989. Biology and systematics of ticks (Acari: Ixodida). Annu. Rev. Ecol. Syst. 20, 397Â…430. Perkins, S.E., Cattadori, I.M., Tagliapietra, V., Rizzoli, A.P., Hudson, P.J., 2006. Localized deer absence leads to tick ampliÂ“cation. Ecology 87, 1981Â…1986. Prvot, C., Licoppe, A., 2013. Comparing red deer ( Cervus elaphus L.) and wild boar ( Sus scrofa L.) dispersal patterns in southern Belgium. Eur. J. Wildl. Res., 1Â…9. Putman, R.J., 1986. Foraging by roe deer in agricultural areas and impact on arable crops. J. Appl. Ecol 23, 91Â…99. Putman, R.J., 1988. The Natural History of Deer. Christopher Helm, London. Putman, R.J., Moore, N.P., 1998. Impact of deer in lowland Britain on agriculture, forestry and conservation habitats. Mamm. Rev. 28, 141Â…164. Randolph, S.E., 2008. Tick-borne encephalitis incidence in Central and Eastern Europe: consequences of political transition. Microbes Infect. 10, 209Â…216. Randolph, S.E., Craine, N.G., 1995. General framework for comparative quantitative studies on transmission of tick-borne diseases using Lyme borreliosis in Europe as an example. J. Med. Entomol. 32, 765Â…777. Reimoser, F., Putman, R., 2011. Impacts of wild ungulates on vegetation: costs and beneÂ“ts. In: Putman, R., Apollonio, M., Andersen, R. (Eds.), Ungulate Management in Europe: Problems and Practices. Cambridge University Press, Cambridge, UK/New York, pp. 144Â…191. Rizzoli, A., Hauffe, H.C., Tagliapietra, V., Neteler, M., Ros, R., 2009. Forest structure and roe deer abundance predict tick-borne encephalitis risk in Italy. PLoS ONE 4, e4336. Ros, R., Pugliese, A., 2007. Effects of tick population dynamics and host densities on the persistence of tick-borne infections. Math. Biosci. 208, 216Â…240. Ruiz-Fons, F., Gilbert, L., 2010. The role of deer as vehicles to move ticks, Ixodes ricinus between contrasting habitats. Int. J. Parasitol. 40, 1013Â…1020. Skalski, G.T., Gilliam, J.F., 2003. A diffusion-based theory of organism dispersal in heterogeneous populations. Am. Nat. 161, 441Â…458. Smit, R., Bokdam, J., den Ouden, J., Olff, H., Schot-Opschoor, H., Schrijvers, M., 2001. Effects of introduction and exclusion of large herbivores on small rodent communities. Plant Ecol. 155, 119Â…127. Sunde, P., Olesen, C.R., Madsen, T.L., Haugaard, L., 2009. Behavioural responses of GPS-collared female red deer Cervus elaphus to driven hunts. Wildl. Biol. 15, 454Â…460. Tagliapietra, V., Ros, R., Arnoldi, D., Cagnacci, F., Capelli, G., Montarsi, F., Hauffe, H.C., Rizzoli, A., 2011. Saturation deÂ“cit and deer density affect questing activity and local abundance of Ixodes ricinus (Acari, Ixodidae) in Italy. Vet. Parasitol. 183, 114Â…124. Talleklint, L., Jaenson, T.G.T., 1997. Infestation of mammals by Ixodes ricinus ticks (Acari: Ixodidae) in south-central Sweden. Exp. Appl. Acarol. 21, 755Â…771. Trdan, S., Vidrih, M., 2008. Quantifying the damage of red deer ( Cervus elaphus ) grazing on grassland production in southeastern Slovenia. Eur. J. Wildl. Res. 54, 138Â…141.
S. Li et al. / Ecological Modelling 276 (2014) 1 Â… 13 13 Vernesi, C., Pecchioli, E., Caramelli, D., Tiedemann, R., Randi, E., Bertorelle, G., 2002. The genetic structure of natural and reintroduced roe deer ( Capreolus capreolus ) populations in the Alps and central Italy, with reference to the mitochondrial DNA phylogeography of Europe. Mol. Ecol. 11, 1285Â…1297. Vor, T., Kiffner, C., Hagedorn, P., Niedrig, M., Ruhe, F., 2010. Tick burden on European roe deer ( Capreolus capreolus ). Exp. Appl. Acarol. 51, 405Â…417. Wahlstrom, L.K., Liberg, O., 1995. Patterns of dispersal and seasonal migration in roe deer ( Capreolus capreolus ). J. Zool. 235, 455Â…467. Wang, H.-H., Grant, W.E., Teel, P.D., 2012. Simulation of climateÂ…hostÂ…parasiteÂ…landscape interactions: a spatially explicit model for ticks (Acari: Ixodidae). Ecol. Model. 243, 42Â…62. Watts, E.J., Palmer, S.C.F., Bowman, A.S., Irvine, R.J., Smith, A., Travis, J.M.J., 2009. The effect of host movement on viral transmission dynamics in a vector-borne disease system. Parasitology 136, 1221Â…1234. Zhang, X.X., Johnson, S.N., Crawford, J.W., Gregory, P.J., Young, I.M., 2007. A general random walk model for the leptokurtic distribution of organism movement: theory and application. Ecol. Model. 200, 79Â…88.