Modelling the effects of habitat and hosts on tick invasions

AbstractMany tick species are invading new areas because of anthropogenic changes in the landscape, shifting climatic variables and increasing populations of suitable host species and tick habitat. However, the relative influences of habitat and hosts in tick dispersal and tick population establishment remain in question. A spatially explicit agent-based model was developed to explore the spatio-temporal dynamics of a generic tick population in the years immediately following the introduction of ticks into a novel environment. The general model was then adapted to investigate a case study of two recent tick species invasions into the Mid-Atlantic United States. The recent simultaneous range expansions of two ixodid tick species, Ixodes affinis and Amblyomma maculatum, provided an opportunity to determine if invasion patterns observed in the field could be replicated in silico on a small scale. The models presented here indicated that for generalist parasites, habitat connectivity is a better indicator tha...


Introduction
Ticks are blood-feeding ectoparasites that parasitize humans and animals and are second only to mosquitoes in spreading vector-borne diseases worldwide (Dennis, Goodman, & Sonenshine, 2005). Many tick species are expanding their ranges as a result of anthropogenic changes in the landscape, shifting climatic variables and increasing populations of suitable host species and suitable tick habitat (Childs & Paddock, 2003;Ogden et al., 2008Ogden et al., , 2008. Climate change has been forecasted to lead to an overall increase in tick habitat in the coming years and is already facilitating tick range expansions worldwide, leading to increasing disease risks (Cumming & VanVuuren, 2006;George, 2008;Leger, Vourch, Vial, Chevillon, & Mccoy, 2012). It is essential to understand factors that limit tick distribution in order to predict disease emergence, as eradication of ticks and their associated pathogens, once established, may be impossible (Cumming & VanVuuren, 2006;Leger et al., 2012).
Tick invasions are different from invasions by other taxa because tick life history is sharply demarcated between periods of movement on-host and longer relatively stationary periods off-host. In order to understand the movement patterns of ticks across a landscape, the suitability of both abiotic and biotic factors must be considered. Ticks have a complex life history that may involve differing host preferences throughout ontogeny and must find suitable hosts at each life stage in order to feed, grow and reproduce. Ticks depend on the large-scale movements of their hosts to transport them across a landscape and are particularly vulnerable to environmental pressures, such as desiccation, when they are freeliving off-host (Leger et al., 2012). Host specificity is key to any parasites' ability to disperse across a landscape and invade new areas (Kruse, Hare, & Hines, 2011). Because many human-biting species of ticks are generalists and can feed on a variety of avian, mammalian and reptilian species throughout ontogeny, these tick range expansions are likely limited predominantly by environmental and climatic variables, including landscape use, habitat availability and the presence of suitable micro-climates (Cumming & VanVuuren, 2006). Ticks are strongly dependent on both host availability and environmental factors for their survival and reproduction in any habitat (Leger et al., 2012), but the relative importance of hosts and habitats in tick range expansions has never been fully explored.
Models have been used to elucidate the complex life history of ticks and to mitigate tick-borne disease risk. Differential equation-based, age-structured difference and matrixbased models have provided insight into the population dynamics of ticks and the dynamics of tick-borne disease (Gaff, Gross, & Schaefer, 2009;Haile & Mount, 1987;Mount, Haile, & Daniels, 1997;Ros & Pugliese, 2007;Sandberg, Awerbuch, & Spielman, 1992). Spatially explicit components have been added using remote sensing, GIS and partial differential equation models (Bunnell, Campbell, & Squires, 2004;Diuk-Wasser et al., 2010;Radcliffe & Rass, 1984). While helpful, most of these models focus primarily on proportional interactions between ticks and hosts that inform our understanding of tick populations and pathogens, but not individual movement. Spatially explicit agent-based models simulate the actions of individual ticks and hosts and can be used to capture the mechanistic phenomena underlying individual episodes of range expansion (Gaff, 2011;Gaff & Nadolny, 2013;Madhav, Brownstein, Tsao, & Fish, 2004;Wang, Grant, & Teel, 2012;Wang et al., 2015).
Here, a spatially explicit agent-based model was developed to simulate the spatiotemporal dynamics of three-host tick populations in the years immediately following the introduction of ticks to a novel environment. Using this model, derived from the TICKSIM model (Gaff, 2011;Gaff & Nadolny, 2013), it was possible to determine the relative strength of influence that host and habitat-based parameters have on invasion rate, population density, geographic pattern of tick invasion and the genetic diversity in resulting tick populations. In addition to addressing broad questions, the model was used to investigate the case study of two recent tick species invasions in Virginia. The recent simultaneous range expansions of the two ixodid tick species, Ixodes affinis and Amblyomma maculatum, into the Mid-Atlantic region of the US provided an opportunity to compare the relative influences of host and habitat choice on invasion dynamics and genetic connectivity Nadolny et al., 2015;Wright et al., 2011). Information gleaned from literature values, field studies on tick ecology and lab studies on tick genetic connectivity was used to parameterize the model, and the emergent properties were compared to determine if the invasion patterns seen in the field could be replicated in silico on a small scale.
In the sections that follow, we will provide some background information and modelling considerations and describe the model following the protocol recommended for agentbased models by Grimm et al. (2010). The general performance of the model will be evaluated, as will the sensitivity of simulated tick dynamics to changes in habitat and hostrelated parameters. Finally, the applications of the model will be demonstrated through determining the relative influence of habitat suitability and host density on simulated invasions by a generalized tick, and the influence of habitat connectivity on the genetic signatures of newly established populations of simulated I. affinis and A. maculatum.

Background information and modelling considerations
Two tick species, I. affinis and A. maculatum, are concurrently expanding their ranges into the Mid-Atlantic region of the US and have been observed invading in different geographic patterns, and with different genetic signatures (Nadolny et al., 2015).Ixodes affinis has been implicated in the sylvatic cycle of Borrelia burgdorferi, the agent of Lyme disease, while A. maculatum is a known vector of numerous pathogens of medical and veterinary importance, including Rickettsia parkeri, the agent of Tidewater spotted fever (Oliver, 1996;Teel, Ketchum, Mock, Wright, & Strey, 2010). Ixodes affinis is generally found in disturbed forested habitat and is a generalist tick species that feeds on small mammals and birds during immature life stages, and medium and large mammals during the adult stage. This tick species exhibits genetically well-mixed populations that are likely created and maintained through short-distance dispersal events throughout the contiguous forested habitat that is abundant in the Mid-Atlantic (Nadolny et al., 2015).
Amblyomma maculatum is another generalist tick species that feeds on birds and mammals but is found only in disturbed open habitats, which are patchily distributed throughout the Mid-Atlantic (Harrison et al., 2010;Wright et al., 2011). Populations of this tick species are genetically isolated from other nearby populations, and each population is likely founded by multiple long-distance founding events and then maintained by the high densities of rodent hosts that are present in grass-dominated habitats (Nadolny et al., 2015). Both I. affinis and A. maculatum generally complete their life cycle in one year and have significant overlap in the hosts parasitized at all life stages (Harrison et al., 2010;Teel et al., 2010). One notable difference is that I. affinis are not known to feed on domesticated artiodactyls, such as cattle, whereas A. maculatum will readily feed on cattle as adults; because the cattle industry is far less developed in the Mid-Atlantic than in other areas where these tick species are established, cattle and cattle pasture are not explicitly included in our models. Their different range expansion patterns and differences in genetic connectivity can likely be explained by the disparate habitat needs of these tick species rather than differences in host preferences. We hypothesize that differential survival in different habitat types is an important factor in determining genetic and spatial spread of ticks, and that these patterns observed in situ can be modelled through the inclusion of heterogeneous habitats and species-specific mortality rates associated with different habitats.
While genetic connectivity has long been used as a proxy for measuring species dispersal among habitat patches (Ibrahim, Nichols, & Hewitt, 1996), there has been no inclusion of genetic parameters in models describing ticks. It has been suggested that incorporating population genetics into agent-based models would be useful for describing many evolutionary processes (Deangelis & Mooij, 2005), but inclusion of genetic components in agentbased models of range expansion and invasion of any species is rare (Bialozyt, Ziegenhagen, & Petit, 2006;Kekkonen, Wikstrm, & Brommer, 2012;Pertoldi & Topping, 2004). By including genetics in agent-based models of species undergoing range expansions, it is possible to validate models using the genetic diversity and connectivity observed in the field.
The agent-based model described here is derived from previous TICKSIM models, which modelled tick-host interactions and emergent patterns of disease prevalence (Gaff, 2011;Gaff & Nadolny, 2013). The current model has been altered in some significant ways from these previous iterations. First, the presence of a pathogen passed between ticks and hosts has been removed in order to generalize the model beyond a specific tickpathogen system and to focus specifically on tick range expansions. Pathogen dynamics can be reintroduced in later, more complex models. Second, while the initial TICKSIM only tracked ticks and hosts, our version of the model also tracks spatially explicit tick populations, and their appearance in space and time, by colour-coding habitat cells based on the presence or absence of ticks. This allows measurement of invasion rate, spatial pattern of invasion and tick population densities overall and in specific habitats. Third, this model includes heterogeneous habitats and includes a desiccation parameter that affects ticks directly to induce mortality in poor-quality habitat patches. Finally, this model includes maternally inherited genetic haplotypes to simulate patterns of mitochondrial gene flow among tick populations.
Other recent agent-based models have included heterogeneous habitats, as well as multiple hosts which ticks can interact. Recent models of lone star tick (Amblyomma americanum) populations in Texas used tick-host-climate-landscape interactions to simulate field conditions and determine the influence of climate change and seasonality on tick populations (Wang et al., 2012(Wang et al., , 2015. These models predicted tick density increases after the addition of a greenbelt to a Texas city and changes in tick densities with the effects of climate change on the seasonal activities of tick hosts. While both these models and the present model include heterogeneous habitat, multiple host types, host home ranges and climate variables; the present model differs in several important ways. The focus of the present model is on tick invasions, not on established populations of ticks. Although the landscape in the present model is markedly less complex than the landscapes modelled by Wang et al. (2012Wang et al. ( , 2015, each individual tick is tracked, each host movement is tracked and there is higher temporal resolution. One other recent model that examines tick range expansions focused on the role of different host types on the range expansion of the blacklegged tick (Ixodes scapularis) into Canada (Madhav et al., 2004). Madhav et al. (2004) found that long-range hosts (e.g. deer) increase invasion rate, high densities of short-distance hosts (e.g. mice) can slow invasions and migratory birds play an important role in the movement of these ticks across landscapes. In this cellular automata model, the authors measured only the effects of varying host parameters on area colonized, using a simplified, spatially explicit landscape. While our model shares many commonalities with the I. scapularis model, Madhav et al. (2004) did not test the influence of different habitat types on invasions, nor did they monitor invasion rate, tick population connectivity or geographic patterns of tick populations. The I. scapularis model was also deterministic, rather than stochastic, which does not take into account individual tick and host interactions, and modelled the influence of tick burden on invasion, which was not tested in the model presented here. A final difference is that the I. scapularis model operated at a coarse resolution of 1 km 2 , while our model operated at a fine resolution, with a total extent of only 165 ha.
This model is based on the premise that interactions between individual ticks, their hosts and their habitats generate the patterns observed in tick range expansions. By varying host and habitat parameters and including stochastic effects, it is possible to determine the relative influence of host and habitat parameters on tick invasions at the local scale. By measuring invasion rate, tick population density, geographic patterns of tick population establishment, and genetic diversity and connectivity of tick populations, it is possible to answer the following questions: (1) Does host density or habitat quality have the greater influence on tick invasions and (2) How does habitat connectivity influence the genetic connectivity and genetic diversity of invading ticks?

Model description
This model description follows the Overview, Design concept and Details (ODD) protocol for describing agent-based models developed by Grimm et al. (2010) and consists of six elements. The first three elements provide an overview, the fourth element explains general concepts underlying the model design and the last two elements provide details. The following description is for a set of complementary models, a general model (Model S1) and a case study model (Model S2). Model S1 was a general model of tick population establishment, with no tick species-specific inputs, and tick-host interactions simulated within a homogeneous habitat. Model S1 was used to assess general model performance, perform sensitivity analyses and address questions on the influence of host density and habitat quality on tick invasions. Model S2 utilizes the same underlying mechanics as Model S1 (described in detail below), but allows for the investigation of a case study using species-specific inputs from I. affinis and A. maculatum. Model S2 incorporates multiple habitats and species-specific survival rates in each habitat type. This allows for investigation into the effects of habitat connectivity on tick invasion patterns and the resulting genetic diversity of new tick populations.
(1) Purpose The purpose of this model is to simulate the spatio-temporal dynamics and genetic diversity of new tick populations after an initial introduction event to a novel area in response to varying host and habitat parameters, and to better understand the underlying mechanisms leading to the establishment and dispersal of tick populations across a landscape. The results of these simulations will help determine the relative importance of host density, host dispersal distance and habitat suitability in shaping the spatial patterns, invasion rate, population density and genetic diversity and connectivity of newly establishing tick populations. The ability to reproduce the spatial and genetic connectivity patterns observed in situ of invading tick species I. affinis and A. maculatum is of particular interest.
(2) Entities, state variables, and scales (a) Agents/individuals This model considered the interactions among three populations of agents: long-distance dispersing hosts (e.g. deer), short-distance dispersing hosts (e.g. mice) and ticks. Hosts were characterized by the following state variables: identification number, home base, home range size, mortality rate, number of ticks currently feeding on the host and the maximum number of ticks able to attach to the host at one time. To keep host populations constant, if a host died it was immediately replaced by another host, which was created on a random cell. The home base of each host was the X, Y coordinates of the cell it was created on. Each host had a specific home range and was only able to move within a certain subset of cells away from their home base. Host categories varied in the distance they could travel per time step and the size of their home range. Hosts moved ticks that were attached across the landscape, and ticks could only move when on a host. A host could carry up to the specified maximum number of ticks, and if a host died, all ticks on that host also died. Ticks were characterized by the following state variables: identification number, sex, life stage, activity, identity number of current host and maternally inherited genetic haplotype. Ticks were assigned a sex (male or female) at birth and moved through the following four life stages throughout ontogeny: egg, larva, nymph and adult. Tick host preferences changed depending on life stage and were reflected by probabilities of successful attachment to each host category. Ticks moved through three activities during each life stage: resting (which includes developing), questing and feeding. Adult female ticks completed a final activity, laying eggs, after feeding. The tick population did not remain constant. Mating was not explicitly included in this model, but ticks were assumed to mate onhost, so female ticks were able to lay eggs after a successful bloodmeal. There was a set number of haplotypes divided equally between the initial ticks at the start of each simulation (e.g. if there were eight initial haplotypes and 32 initial ticks, there would be four ticks of each haplotype), and each new tick 'hatched' throughout the course of the simulation inherited its haplotype from its mother. Parameter values for hosts and ticks can be found in Table 1. Parameter values were derived from the literature where available, or parameter values were estimated for a generalized tick model (GTM) based on literature reviews and field data. (b) Spatial units Environmental conditions did not change on a cell by cell basis, but cell colour was used to indicate changes in tick occupancy patterns. Cells were either unoccupied, occupied (between one and five ticks of any life stage present) or populated (six or more ticks present on that cell). New tick populations were observed through the colour variables of cells, with colour reflecting patch occupation by ticks. A green colour variable indicated the background environment where no ticks are present. Once a tick was hatched or moved onto a cell, that colour variable turned to yellow to indicate that ticks were present. If six or more adult ticks occupied a cell simultaneously, that colour variable changed to red to indicate that a population of ticks was present  Notes: LD indicates long-distance dispersing hosts (e.g. deer), SD indicates short-distance dispersing hosts (e.g. mice). A reference or a reason for the assumption is provided for each parameter. SA indicates that values were chosen based on computational limits determined by sensitivity analyses or trial and error. GTM indicates that values were chosen to represent a generic three-host tick model, with larvae and nymphs feeding primarily on small mammals and adults feeding on deer, informed by lists of preferred hosts for I. affinis (Harrison et al., 2010) and A. maculatum (Teel et al., 2010). (Fish & Howard, 1999). Once the colour variable for a cell had changed from green to yellow or from yellow to red, it would remain changed until the first day of the next year, when all cells were reset to green. The occupancy of that cell did not change, only the colour variable; any cell that was either occupied or populated would immediately turn back to the appropriate yellow or red colour during the first time step of the new year. This enabled output at the end of the simulation to reflect only the most recent years tick occupancy patterns across the simulated landscape. We chose to reflect only the most recent years tick occupancy patterns because when sampling a tick population in situ, only the ticks present are able to be sampled. Cell colours were recorded at each time step, to be used as a proxy for tracking tick population establishment. Including multiple years worth of ticks in our final sampling output would both inflate the number of ticks and unreasonably expand their area of occupancy. The presence or absence of ticks of each haplotype across space was also monitored by a cell variable. The haplotype variable recorded if there were any ticks of each haplotype as a simple presence/absence variable for each cell. Like the colour variable, once a haplotype had been recorded in a cell, that record remained until the first day of the next year, when all haplotype variables were reset to enable only the most recent years distribution of haplotypes across the landscape to be recorded at the end of the simulation, for reasons described above. (c) Environment The environment was set up as a grid of 51 × 51 patches, with each cell roughly representing .06 ha for a total simulated area of roughly 165 ha, with hard (reflective) boundaries. Hard boundaries were chosen because wrapping boundaries would result in unrealistic jumps in tick occupancy from one end of the simulation to the other. The highest hierarchical level in the model was the abiotic environment and its fluctuations. Type of habitat was determined by a desiccation parameter D that influenced tick survival at all life stages when the ticks were off-host (questing or resting). Desiccation could be increased to increase tick mortality (indicating habitat of poorer quality), or decreased to not influence tick mortality (indicating good habitat where ticks were easily able to survive off-host). The general model was run with one, homogenous habitat (one desiccation parameter for the whole environment, Model S1), but we incorporated multiple habitats in different parts of the environment, each with its own desiccation parameter to investigate the effects of habitat connectivity on tick invasions (Model S2). Regardless of whether a single habitat or multiple habitats with different desiccation values were used, habitats were constant and the desiccation parameters did not change during the simulation.
Time of year (t) also factored into tick mortality. Each time step represented one day, and ticks were more likely to die in late fall (October and November), winter (December through March) and mid-summer (July) than in other months when weather conditions were more favorable (Table 1). Each simulation ran for 1080 time steps, or 3 years, to give tick populations sufficient time to establish and spread.
(3) Process overview and scheduling The model proceeds in daily time steps. Within each day or time step, six modules happen in the following order: set day of year, tick changes, host changes, calculate new populations, calculate haplotype distributions and reset ( Figure 1). Within each module, individuals are processed in random order. (4) Design concepts (a) Basic principles The underlying principle of this model is that independent agents interact with one another and simulate the interactions that ticks would have with hosts in the field. Ticks interacted with hosts using them as blood meals, and ticks had a probability of attaching to each host type that changed throughout ontogeny, with higher probabilities indicating increased preference for that host. Each tick was assigned a genetic haplotype that was passed onto offspring. Hosts moved across the landscape according to their type (i.e. species) and were constrained by home ranges. Tick mortality was influenced by habitat (desiccation), time of year and host availability. Through these interactions, ticks were transported across the environment by hosts and established new populations with specific genetic signatures.
The emergent property being modelled is the establishment of new populations of ticks that then either sustain themselves or die out. Population establishment was measured in four ways: invasion rate, tick population density, genetic diversity and connectivity, and patterns of spatial spread. (c) Sensing Ticks sensed hosts only within their own cell and had a given probability of successful attachment and feeding on that host (Table 1). After host movement in each time step, ticks sensed hosts and either attached successfully or did not successfully attach and continued to quest. Ticks could not attach to hosts that were already carrying the maximum number of ticks for that host. Ticks also interacted with the environment, as both desiccation and time of year were factored into calculating tick mortality during each time step. Hosts did not sense or interact with other hosts and moved around the environment independently of one another in this simple model. Ticks did not sense or interact with other ticks, even to mate; any adult female that had successfully completed a bloodmeal would reproduce, as ticks mate on-host and probability of predation on fed females was not included in this model.

(d) Interaction
Ticks sensed hosts within their cell and attached to a host to attempt to successfully obtain a blood meal. Once a tick had successfully attached to a host, it switched from 'questing' activity to 'feeding'. A tick would feed on a host for a specific number of days. Once the tick was engorged, it would detach from the host and drop off on whatever cell that host had moved to during the interim time steps. The tick would then transition to 'resting', while it transitioned from one life stage to the next. After a given number of time steps, the tick would again begin questing and would attempt to attach to its next host. Adult female ticks laid a set number of eggs after completing their final bloodmeal. Through the interactions of ticks with hosts, ticks could be transported across the landscape and seed new populations with specific genetic signatures.
(e) Stochasticity Stochasticity was included in calculating host movement, host mortality, tick mortality and successful tick attachment to hosts while questing. Probabilities for stochastic parameters are described in Table 1. The processes were stochastic for each run and each agent had equivalent fitness. (f) Observation The following metrics were monitored throughout each time step: the number of ticks in each life stage, the number of occupied (yellow) and populated (red) cells, the location of each occupied cell, the location of each populated cell, the number of ticks of each haplotype and the presence or absence of ticks of each haplotype in each cell. Cell colours were used as a proxy for tracking tick population establishment or tick presence throughout year. Metrics were used to determine the following four outputs: tick population density, tick invasion distance, tick genetic diversity and the spatial pattern of the tick spread. Tick population density was measured by two metrics, the number of cells occupied by at least one tick (yellow cells) and the number of cells with tick populations (six or more ticks present, red cells). Tick invasion distance was measured using four outputs: the distance from origin to the furthest cell occupied by at least one tick, the distance from origin to the furthest cell containing a population of ticks, the mean distance from the origin to occupied cells and the mean distance from the origin to cells containing a tick population. Tick genetic diversity was measured by the number of surviving haplotypes. Finally, the spatial pattern of tick spread was measured by the number of clusters of occupied and populated cells. A cluster was defined as being three or more cells away from the nearest red or yellow cell, and made of two or more red or yellow coloured cells.

(5) Initialization
Each simulation began with a landscape with a set number of randomly distributed hosts of both species, and an initial number of ticks created in the cell at the centre of the simulation (at X,Y coordinate 0,0). Each simulation ran for 1080 time steps, or 3 years, to give tick populations sufficient time to establish and spread. The initial state of the simulation used base parameters outlined in Table 1, which were derived from values in the literature, and in previous versions of this model (Gaff, 2011;Gaff & Nadolny, 2013 (Mount et al., 1997;Randolph, 1994).
(b) Process tick life cycle Each time step, each individual tick had a given probability of dying. Tick mortality (M) depended on time of year (t), habitat (represented by the desiccation parameter D, which could change depending on the type of habitat ticks were in, as in Model S2), total tick population size (N) and carrying capacity, or maximum number of ticks permitted in the simulation (k), and was calculated as where t reflected the relevant value from Table 1; N was re-counted at the beginning of each time step; D was a constant that could be varied from 1 (no effect on tick mortality) to 50 (significant influence on tick mortality) and k was set to 50,000 as a computational constraint and a reasonable environmental limitation. Once per time step, a random number between 0 and 1 was generated for each tick; if the random number was less than M, the tick would die. Ticks could rest (indicating time spent molting or in quiescence between life stages), quest, feed and reproduce. After a certain amount of time in a given activity, ticks would change to the next chronological activity (Figure 1).Ticks quest to find hosts, take bloodmeals and progress to the next life stage. If a tick exceeded the maximum questing time (Table 1), it would die. Each tick in a questing life stage (larvae, nymphs and adults) began questing by identifying a potential host, which was any host occupying the same habitat cell with fewer than the number of maximum ticks per host already attached. If a potential host was present, ticks would successfully attach if a random number generated between 0 and 1 was less than the attachment probability specific to that host and life stage (Table 1). For example, larval ticks had a higher likelihood of successfully attaching to a mouse than a deer. After successfully attaching, lists were updated to track the number of ticks currently feeding on each host, and the ticks would feed for six days. After feeding, the fed tick would drop off in whatever habitat cell the host was in at the end of the six days and would change activity to the next activity laid out in Figure 1. Adult males would die after feeding and presumably mating, while adult females would lay eggs and then die. It is assumed that adult males of both species fed only once, on one hostif they drop off, they die. (c) Process host mortality and movement Long-distance and short-distance host movement and mortality was processed the same way, albeit with different numbers allocated to distance and home range size (Table 1). When hosts were created, the centre of their home range would be the habitat cell in which they were 'born'. Each time step, hosts could move a given number of cells (deer or mouse rate, Table 1) in any direction from the centre of their home range. If the host was farther than a given number of cells (deer or mouse home range, Table 1) away from the centre of their home range, they would turn to move back towards that centre during the next time step. In this way, hosts would always remain within a home range of a given size. Each time step, each individual host had a given probability of dying. Host mortality was set at .02 (Table 1) for both short-distance and long-distance hosts to maintain some level of turnover in the simulation. Each time step, a random number between 0 and 1 was generated; if the random number was less than .02, that host would die. To keep host populations constant, for each host that died another was created immediately and placed on a random habitat cell. If a host died, all ticks on that host also died, but were not replaced.

(d) New populations
In homogenous habitat (Model S1), populations of ticks were visualized by changing cell colour. To visualize only recent tick populations, the visualization of populations would reset after every 360 time steps, as described above. If there were no ticks on a habitat cell, the cell would remain its green background colour. If there were one or more ticks on a cell, the cell would change colour to yellow to indicate the cell was occupied. If there were six or more adult ticks on a cell, the cell would turn red, to reflect the presence of a population, i.e. that the cell was populated. From then on, the number of ticks on each cell would be counted every time step, and cells would change colour to reflect tick occupancy. Cells would stay coloured until the next reset at the beginning of the next year.
In heterogeneous habitat (Model S2), populations were not visualized, but the number of ticks present in each cell were recorded. (e) Determine genetic connectivity To keep track of genetic connectivity, the number of ticks of each haplotype was counted on each cell and updated each time step. This resulted in a running total of the maximum number of ticks with each haplotype that had been present on each cell. This number was reset each year, so that only the most recent occupancy since the beginning of the year was recorded for each cell.
The model was programmed using NetLogo version 5.0. This software was written by Uri Wilensky in 1999 and is freely available (http://ccl.northwestern.edu/netlogo/). The model was run on the high-performance computing cluster using BehaviorSpace to run experiments 'headless', from the command line. All statistical analyses on model results were completed in R and MATLAB (MathWorks, Inc, 2015;R Core Team, 2015).

General model performance
To evaluate general model performance, the model was initialized with 32 nymphal ticks of eight haplotypes, four ticks of each haplotype, dropped at the centre of homogenous habitat with a desiccation parameter of 1.0, which does not influence tick mortality. Tick mortality was therefore influenced only by seasonal and stochastic effects. The simulation was started on 01 March, to begin in the spring when nymphs of I. affinis and A. maculatum are active. All parameter values were set at base values as indicated in Table 1. The resulting phenologies of each tick life stage (Figure 2, based on an average of 25 runs with Table  1 inputs) show nymphs peaking in the spring, adults peaking in the summer and larvae peaking in the winter each year. Growing annual numbers of nymphs and adults indicated that the tick population was increasing every year, and that the model was coded properly.

Sensitivity analyses
To identify the factors with the greatest influence on tick invasion, sensitivity analyses were performed on nine parameters ( Table 2). Each of the parameters was varied one at a time as shown in Table 2, while keeping all other parameters constant as shown in Table 1. Each variant of each parameter was run 25 times until 1080 time steps (3 years), or until the tick population was extirpated. A priori power analyses conducted with G*Power version 3.0.10 indicated that 25 runs per parameter set would be sufficient to provide statistical significance (Faul, Erdfelder, Lang, & Buchner, 2007). For each parameter variation, metrics to determine the following four outputs were recorded at the last time step of the simulation: tick population density, tick invasion distance, tick genetic diversity and the spatial pattern of the tick spread. Cluster analysis was computed using the DBSCAN package for MATLAB (MathWorks, Inc, 2015). Spearman rank correlation coefficients were calculated to determine the relative influence each parameter exerted on the measured variables (Tables 3 and 4). Some patterns that emerged from the sensitivity analyses were striking, though not wholly unexpected. Mouse density and mouse home range size exerted the largest influences on tick population density at the end of the third year, with more mice with larger home ranges resulting in higher tick densities. Deer home range size had the greatest influence on invasion rate and on the spatial patterns of tick invasions in heterogeneous habitat. Larger home ranges of long distance dispersers resulted in faster spread and more clustering. The maximum number of ticks attached to a deer was not significantly correlated with any invasion metric, although increasing the number of ticks that could be attached to a mouse positively influenced the number of populations established. The initial number of haplotypes had the greatest influence on the percentage of haplotypes surviving at the end of three years. More initial haplotypes resulted in more die-off of haplotypes through stochastic effects and a smaller percentage surviving through the end of the third year. However, increasing the number of initial haplotypes significantly increased the mean number of haplotypes surviving at the end of year three (Spearman correlation coefficient = .96, p < .01).

General model application: hosts density vs. habitat quality
One essential question when dealing with tick range expansion is the relative influence of habitat and hosts in determining dispersal and establishment. Because ticks spend much of their lives off-host, generalist ticks may establish anywhere that has suitable climatic conditions (Cumming & VanVuuren, 2006). However, ticks depend on hosts not only for sources of bloodmeals, but also to move them across a landscape.
To test whether host density or habitat quality had the largest influence on tick invasions, a set of four models was run 100 times per parameter set. The models varied habitat quality between good and poor, and host density between high and low (Figure 3). Changing the number of available mice and deer varied host density, and changing the desiccation parameter varied habitat quality. High host levels were set at 1600 mice and 120 deer, and low host levels were set at 200 mice and 15 deer. Good quality habitat was created by a low desiccation parameter that did not influence tick mortality (D = 1), while introducing a high desiccation parameter (D = 30) created poor quality habitat. The magnitude difference between high and low desiccation parameters was larger than that between high and low host densities, because humidity can vary tremendously between adjacent habitat fragments, but adjacent host densities exhibit lower variation because of host movements. Each model was run on a landscape of homogenous habitat that did not change over time, and all metrics other than host density and habitat were kept at baseline levels (Table 1).  Notes: Trends depicted are averages from 25 runs of the general tick Model S1 using the parameters outlined in Table 1.
Invasion metrics were measured at the last time step after 3 years for each run. Partial rank correlation coefficients (PRCC) were computed to determine the relative influence of hosts and habitats on invasion metrics (Table 5), and boxplots were generated to display the differences between means among the four models (Figure 4). In homogenous habitat, host density was found to have a relatively greater influence on almost all invasion metrics, although both host density and habitat quality showed significant effects on all invasion metrics (Table 5). Higher host densities resulted in higher tick densities, faster invasion rate, more clusters and more genetic diversity of ticks at the Notes: Habitat quality was modelled as either good? (D = 1) or bad? (D = 30), and host density was either high? (1600 mice, 120 deer) or low? (200 mice, 15 deer). These images represent example runs of the tick populations exposed to these conditions across homogeneous landscapes after three years. Green cells are unoccupied habitat, red cells indicate populated cells with > 6 adult ticks, and yellow cells are occupied by at least one tick of any life stage. end of three years. Higher desiccation, or poor habitat, resulted in lower tick density, slower invasion rate, fewer population clusters and less genetic diversity. Interactions between host quantity and habitat quality produced different outcomes across invasion metrics ( Figure  4). For example, when hosts were abundant and habitat quality was high, invasion distance, density and genetic diversity were all higher than in any of the other models. The number of clusters, however, was highest when hosts were abundant but habitat was poor, indicating that tick populations likely established for only short durations before winking out, but ticks were able to feed, move and reproduce freely enough that many small populations could remain for at least three years. The number of tick population clusters was lowest when populations were only rarely able to persist (few hosts and poor habitat), but the number of clusters was also low when one large population resulted in one large cluster across the whole simulated environment.

Case study model application: influence of habitat connectivity on genetics (Model S2)
While the above results suggested that host density is more important than habitat quality, field data can show otherwise. Since 2009, our lab has conducted an intensive tick field surveillance effort in southeastern Virginia, during which we observed two tick species invading across the landscape (Nadolny, 2016;Nadolny, Wright, Sonenshine, Hynes, & Gaff, 2014;Nadolny et al., 2011;Wright et al., 2011). We found that although I. affinis and A. maculatum were expanding their ranges simultaneously, they were doing so in different spatial patterns and with different genetic connectivity and structure in the resulting populations (Nadolny, 2016;Nadolny et al., 2015). We found I. affinis in contiguous populations in the well-connected disturbed forested habitats that are common in the southeast, and I. affinis populations exhibited high genetic homogeneity as a result. Conversely, A. maculatum populations were only found in grassy habitat patches and were often only present for a period of a few years before the process of ecological succession extirpated the local population (Nadolny, 2016). Amblyomma maculatum tick populations found in these habitat patches showed genetic dissimilarity to other populations, with each population being genetically distinct from its neighbours (Nadolny et al., 2015). Despite these differences in spatial and genetic patterns, both tick species are trophic generalists, parasitizing a diverse and overlapping assortment of birds and small mammals as immatures, and feeding on large mammals such as white-tailed deer as adults (Harrison et al., 2010;Heller et al., 2016;Teel et al., 2010).
To address the importance of habitat heterogeneity and test the effects of habitat connectivity on tick genetic connectivity, two models, each with multiple habitat types, were adapted from the generic tick model by varying the desiccation parameter in designated regions of the simulation ( Figure 5). Because both tick species are trophic generalists but appear to be habitat specialists, the same life history parameters (Table 1) were used to model each tick species. The two models were each run 100 times, with one modelling ticks moving through well-connected matrix of high quality habitat with patches of poor habitat and one modelling ticks moving through a matrix of poor habitat to patches of high-quality habitat. The model with well-connected high-quality habitat represented I. affinis moving through contiguous preferred forested habitats and experiencing higher mortality in field patches, while the other model represented the observed tendency of A. maculatum to establish only in patchy field habitats and experience increased mortality in well-connected forested habitats ( Figure 5). The configurations of habitat were kept consistent between the two models, but the desiccation parameters were reversed to indicate each ticks respective habitat preferences. The presence or absence of each haplotype was recorded in each cell, and each habitat type was separately tracked to determine the genetic makeup of the tick population within that habitat at the end of three years. Tracking the presence of haplotypes in each habitat also provided a proxy for tracking tick presence in each habitat. To recreate the patterns of genetic and spatial connectivity, we observed in the field and we would expect the I. affinis model to result in one contiguous tick population with a similar genetic makeup throughout the forested area of simulation, while the A. maculatum model should result in genetically distinct populations in each patch of field habitat ( Figure 5).
To analyse the similarity between tick populations in each habitat, a binary similarity coefficient was used. These are simple similarity measures based on presence-absence data, Table 4. P values for the Spearman rank correlations from sensitivity analyses shown in Table 3. Note: Occupied patches had one or more ticks on a cell, and populated patches had six or more ticks occupying a cell. Figure 4. Boxplots comparing means between four treatments groups testing influence of host density (with host density designated as 'high' or 'low') and habitat quality (with habitat quality designated as 'bad' or 'good') on invasion metrics. Metrics are (A) tick invasion distance, measured by adding together all rate metrics, (B) tick invasion distance, measured by multiplying together all rate metrics, (C) tick population density, measured by adding together all density metrics, (D) tick population density, measured by multiplying together all density metrics, (E) genetic diversity and (F) number of clusters.

P-values
Notes: Rate metrics included maximum and average distances to occupied and populated patches, and density metrics included number of populated and occupied patches. In lieu of presenting all box plots produced for rate and density metrics, only the additive (+) and multiplicative (x) plots are shown to give an overview of the trends seen. Lines through the middle of each box represent the median, and whiskers extend to lower and upper quartiles (25th and 75th percentile). Outliers are represented with clear circles.
usually of species in habitats (Krebs, 2014). Here, haplotypes were used instead of species, and Sorensens index was used to compute a similarity coefficient between tick populations in each pair of habitat types (Krebs, 2014). Because there were three habitats (one forest and two field patches), the results of all three pairwise comparisons were summed to create a general similarity index that was used for statistical analysis. Because many cells did not have haplotypes present, the data were zero inflated and non-normal, so Wilcoxon rank sum tests with continuity corrections were used to compare means between the treatment Notes: PRCCs are between input parameters of the models (host density vs. habitat quality) and the output metrics measuring tick population density, spatial pattern of tick spread, tick genetic diversity and tick invasion distance. All PRCCs except those indicated with * were significantly different, with p < .01. Additive and multiplicative indices were generated through adding together or multiplying together all tick density metrics (number of occupied and populated patches) or tick rate of invasion metrics (distance to furthest occupied and populated patches, and mean distance to occupied and populated patches).
groups. The general similarity index, mean haplotype survival overall and mean haplotype survival in each habitat type were compared between treatment groups ( Figure 6).
To determine which tick species exhibited populations that were more genetically similar, the general similarity index was compared between the I. affinis model and the A. maculatum model. We found that I. affinis populations were more similar among habitats, while populations in each habitat in the A. maculatum models tended to be more genetically distinct. Habitat connectivity also influenced the number of haplotypes that persisted until the end of the simulation, with I. affinis models resulting in more haplotypes surviving than A. maculatum models. In I. affinis models, there was also a greater diversity of haplotypes in forested habitat, while A. maculatum models resulted in proportionally more haplotypes present in field habitats. These results correspond with population genetics patterns described by Nadolny et al. (2015), where I. affinis across the eastern US comprised two large, contiguous, genetically similar populations, while each sampling site with a population of A. maculatum across the eastern US was genetically distinct.

Discussion
No model is a perfect representation of the agents and environment that it represents, and the present model has several limitations that are fully acknowledged. Because of the stochasticity inherent in agent-based models, they generally do not provide proof of a stable system, there is no single solution and results may not be generalizable or easy to validate (Grimm & Railsback, 2005). However, most of these limitations can be overcome by running sufficient numbers of simulations to avoid outliers. This model takes place on a very small spatial scale, which, while useful for modelling local invasion patterns, Figure 5. Adapted model designs to investigate a case study of Ixodes affinis and Amblyomma maculatum, and simplified tick populations depicting spatial and genetic differences in tick populations caused by habitat preferences. (A) Models for testing the influence of habitat connectivity on genetic connectivity, and reproducing patterns exhibited by I. affinis (left) and A. maculatum (right). In the I. affinis model, initial ticks are dropped in the centre of the completely connected 'forest' habitat, which has a desiccation parameter, D, of 1, indicating good tick habitat. The two 'field' patches (red and yellow) both are poor habitat quality with D = 20 for I. affinis ticks. The A. maculatum model has the reverse setup, with initial ticks dropped in forest habitat that poor quality, with nearby field patches of higher quality. In both simulations, ticks were initialized at the 0,0 coordinate, in the centre of the forest habitat, and life history parameters were identical. (B) Graphic depicting simplified populations of I. affinis (left) and A. maculatum (right), derived from population genetics data published in Nadolny et al. (2015). Yellow habitat patches represent field habitats, and green represents forested habitat, as in the model depicted in 5A. Pie charts represent populations of each tick, with different colours representing different haplotypes present in each population. Arrows represent pathways of dispersal between populations. Tick images suggest where ticks of each species are most likely to be found. Because I. affinis survives better in forested habitat, which is well connected in the Mid-Atlantic United States, and in our simulation, we expect similar genetic makeups in each population, with little survival in field habitats. In contrast, because A. maculatum is better suited to patchy field habitats, we expect each habitat patch to contain its own population, with little genetic similarity between populations and little survival in forested habitats. is unsuitable for broader inferences on a regional scale. The number of individual ticks and hosts that could be included limited the scale of the model. It was decided to scale the model down in order to maintain the individual-based nature of both ticks and hosts, rather than inferring tick-host interactions and movements as other models have done (Wang et al., 2012(Wang et al., , 2015. This model used the standard threshold for an established population of more than six ticks (Fish & Howard, 1999), which is appropriate for these initial, small populations. Also, in this simple model, hosts do not interact with habitat or other hosts, and the host population is kept constant, so host demography is not included. Also, the edges of the model simulation are reflective, which would have strong influences after more than a few years of simulated time. For that reason, simulations were kept to a maximum time of three years, which enabled tick populations sufficient time to establish and grow, but not 'outgrow' the small area of the simulation.
The general tick model yielded several interesting inferences regarding the effects of hosts and habitats on tick range expansions. The results of sensitivity analyses indicated that in homogeneous habitat, long-distance dispersing hosts (e.g. deer) had the greatest influences on both invasion rate and the spatial patterns by which ticks invade, while short-distance hosts (e.g. mice) had the greatest influence on tick densities. This supports the results from a spatially explicit cellular automata model of I. scapularis invasion, which indicated that ticks feeding primarily on deer increased invasion rate, while feeding on mice reduced invasion rate (Madhav et al., 2004). However, our results also highlight the importance of small mammal communities during the initial establishment of ticks, and the link between abundant small mammals and the ability of an area to sustain large numbers of ticks, which has never before been explicitly modelled.
Results from the general tick model comparing the relative influences of host density with habitat quality suggested that both significantly influence tick invasion rate, pattern, density and genetic diversity. Interestingly, in homogenous habitat, host density had the largest influence on all invasion parameters measured. However, once heterogeneous habitat was introduced in the adapted model to recreate patterns exhibited by I. affinis and A. maculatum, connectivity between high-quality habitat patches significantly influenced tick invasions. Habitat connectivity especially influenced the genetic connectivity between tick populations, even at a very small spatial scale. Greater habitat connectivity resulted in more genetic similarity between populations and more overall genetic diversity surviving at the end of three years. Isolated habitat patches resulted in more genetically distinct populations, with less overall genetic diversity surviving after three years.
These changes in genetic connectivity and diversity mimic the genetic patterns exhibited by I. affinis and A. maculatum in the field, where habitat preference resulted in significant differences between the two species (Nadolny et al., 2015). In the model, the I. affinis simulation resulted in one genetically heterogeneous population concentrated in the forest, and the A. maculatum simulation resulted in several genetically distinct populations concentrated in the field patches. These results replicate the genetic patterns found in I. affinis and A. maculatum populations across the eastern US (Nadolny et al., 2015). Simply changing the connectivity between preferred habitat patches for each tick species simulated these genetic and spatial differences for the parameters used in these scenarios. Knowledge of the genetic makeup of these populations in situ provides a measure of validation for this model, and the generation of such close patterns suggests that the agent-based model was robust, even at such a small spatial scale. The only alteration made between the models representing I. affinis and A. maculatum was in the spatial structure of the desiccation parameter; therefore, that is the only condition that could have been responsible for the changes in patterns observed. Although there are other ecological differences between the two tick species that could inform invasion patters, many of the important parameters including common hosts parasitized, length of time feeding and lifespan are markedly similar.

Future directions
The models presented comprise an initial investigation of some parameters which may be important to tick range expansions and serve as a proof-of-concept for the utility of agent-based models to address questions of biological invasions. As a result, there are many things that we did not address. Since the focus of these models are on the initial invasion, we did not track the tick population beyond three years post arrival. A tick carrying capacity was also set to accommodate the computational limitations associated with agent-based models tracking the actions of thousands of individual agents.
Biologically, several important factors should be considered for future iterations of this model. Both I. affinis and A. maculatum parasitize birds as immatures (Harrison et al., 2010;Teel et al., 2010), but due to the computational constraints of maintaining an agentbased model that included individual tick-host interactions; it was impossible to model a landscape of sufficient size to include avian movement. Also, other factors modelled as identical between the two species may actually contribute to differences in spatial and genetic patterns, such as mortality rates, durations of feeding, development periods or feeding behaviour of male ticks (Oliver et al., 1987;Teel et al., 2010). Inclusion of genetic information from male ticks would allow for the use of nuclear chromosomal markers in addition to the maternally inherited markers used in this model. While the general tick model parameters included in Table 1 were within the average range of behaviour for both species, differences in these parameters should be considered in future, more sophisticated models. Biological parameters such as the number of days a tick spent in each life stage or feeding on a host were also left static, whereas it would be more true to life to vary those parameters stochastically.
As written, the model presented here is a useful tool for learning about the relative influence of habitat and host parameters on tick range expansions and genetic connectivity. For future incarnations of this model, it will be useful to measure the spatial patterns of spread more directly, including measuring the size of clusters and distance between clusters. Monitoring the persistence of clusters of tick populations would also provide useful information on the stability of tick populations throughout the years immediately following a range expansion. Including more complex host dynamics, such as herding behaviour, territory defence and interaction of hosts with the landscape will also be important for more realistic models of tick-host interactions.
The effects of changing habitat over time are important for tick population establishment and persistence, so for future models it would be interesting to change the desiccation parameter over time to model the influence of succession or deforestation on tick populations. Also, the simulations presented here were at a very small spatial scale, where edge effects may influence invasion rate after only a few years. Increasing the size of the model, perhaps through the creation of a mixed deterministic/stochastic model, would allow for more realistic model of tick invasions over regional spatial scales. A larger spatial scale could include more host categories, such as migratory birds, and allow a more detailed investigation of the influence of different host categories on the ability of tick populations to establish and persist outside their native range.

Funding
Dr Nadolny was supported by a Science, Mathematics and Research for Transformation (SMART) scholarship from the Department of Defense and the American Society for Engineering Education, as well as a dissertation support grant provided by the Jayne Koskinas Ted Giovanis Foundation for Health and Policy, a Medical, Urban, and Veterinary Entomology (MUVE) grant administered by the Entomological Society of America, and a research grant from the Virginia Academy of Sciences. The content is solely the responsibility of the authors and does not necessarily represent the official views of any granting agency.