Modeling the Impact of Climate and Landscape on the Efficacy of White Tailed Deer Vaccination for Cattle Tick Control in Northeastern Mexico

Cattle ticks are distributed worldwide and affect animal health and livestock production. White tailed deer (WTD) sustain and spread cattle tick populations. The aim of this study was to model the efficacy of anti-tick vaccination of WTD to control tick infestations in the absence of cattle vaccination in a territory where both host species coexist and sustain cattle tick populations. Agent-based models that included land cover/landscape properties (patch size, distances to patches) and climatic conditions were built in a GIS environment to simulate WTD vaccine effectiveness under conditions where unvaccinated cattle shared the landscape. Published and validated information on tick life cycle was used to build models describing tick mortality and developmental rates. Data from simulations were applied to a large territory in northeastern Mexico where cattle ticks are endemic and WTD and cattle share substantial portions of the habitat. WTD movements were simulated together with tick population dynamics considering the actual landscape and climatic features. The size of the vegetation patches and the distance between patches were critical for the successful control of tick infestations after WTD vaccination. The presence of well-connected, large vegetation patches proved essential for tick control, since the tick could persist in areas of highly fragmented habitat. The continued application of one yearly vaccination on days 1-70 for three years reduced tick abundance/animal/patch by a factor of 40 and 60 for R. annulatus and R. microplus, respectively when compared to non-vaccinated controls. The study showed that vaccination of WTD alone during three consecutive years could result in the reduction of cattle tick populations in northeastern Mexico. Furthermore, the results of the simulations suggested the possibility of using vaccines to prevent the spread and thus the re-introduction of cattle ticks into tick-free areas.


Purpose
The model was developed to investigate the effect of anti-tick vaccines on wild whitetailed deer (WTD), Odocoileus virginianus, related to host dispersal in northeastern Mexico. Simulations of the life cycle of the cattle ticks, R. annulatus and R. microplus, was based on already published process-driven models and the effect of vaccination on tick reproductive parameters based on actual vaccine trials. The study was based on simulating a population of WTD inhabiting and moving on a 2D landscape, by using and individual based model where movement is determined by the area of each path of vegetation (carrying capacity) and the distance to near patches.
The model was not used to evaluate a realistic vaccination protocol of WTD to reduce tick infestations in the area because technology for mass vaccination of large wildlife populations is still unavailable for WTD. We rather focused on examining how an actual landscape configuration may interfere with the actions towards reduction of ticks by vaccination, an intervention that must by applied in the complete territory.
The study area covers a region located between 96º and 104ºW and 22º and 29ºN in the States of Coahuila, Nuevo León and Tamaulipas in the northeastern part of Mexico (Fig. 4). Previous simulations of the landscape configuration evaluated the main spread routes of R. microplus from Mexico into the USA [3]. The area is thus known to have a special significance for tick control because it supports large WTD populations [1] that has a demonstrated role in the maintenance of cattle tick populations [4]. WTD in this region are believed to be responsible for the spread of tick populations into the USA [5][6].
The objective of the study was to find the optimum strategy for WTD vaccination against both cattle tick species, which depends on how vaccinated WTD move through the landscape and how they share the same habitat with unvaccinated cattle.
If WTD is more abundant on critical patches of the landscape, they will carry more ticks than cattle and therefore a reduction in the abundance of ticks will be observed.
If the configuration of the landscape does not allows for a critical tick load on WTD, the effects of vaccination will be negligible. This strategy depends not only upon habitat configuration but also on the tick seasonality as driven by actual climate features.
Models were first run in simulated landscapes to assess the impact of the different features of the habitat and relative densities of cattle (unvaccinated) and WTD (vaccinated) on the control of tick populations. We then applied the model to an actual area in northeastern Mexico to examine the probable effects of vaccination on the modelled densities for the two most important tick species affecting livestock and WTD in the area, R. microplus and R. annulatus.

States variables and scales
Animal dispersal in a fragmented landscape depends on the complex interaction between landscape structure and animal behaviour. Our model uses a patch structure, in which the patches are the sites from and to where animals can move, and have variables dimensions. The tick colonization patterns strongly depend on abiotic (climate) suitability, the degree to which the patch is isolated and the overall pattern of interconnections among patches of suitable habitat, which explains host movements and frequency [3]. Previous empirical studies [2,[7][8] support the hypothesis that movements of the main tick hosts among patches support a network of tick dispersal. The movements of the WTD across the network of patches depend only upon the carrying capacity of the patch and simple rules addressing the percentage of animals that will move to another patch. Time of the year (in 10 days intervals) was used to calculate the amount of ticks active and questing for hosts. Therefore, values on tick abundance are used to calculate the amount of ticks carried by hosts in their movements.

Process overview and scheduling
A spatially explicit agent model was developed in which patches of habitat were randomly allocated within a matrix of non-habitat. The model allows the generation of patches within a range of sizes and distances among them to examine the impact of the physical features of the network on animal movements. Encounters between questing ticks and hosts are governed by mass-action. A tick-host encounter results in the transition of the tick to the feeding stage, with a certain probability of moulting success after feeding. The performance of the tick population (e.g. the production of new individuals as output of the current generation) is governed by mechanisms of density-dependent regulation, according to host resistance and acting on tick feeding success. Temperature and relative humidity act on tick survival while moulting or questing on the vegetation. Hosts were assumed to move across the landscape according to the shape and size of every patch and a set of rules involving habitat perception [9]. The density of ticks is regulated by the climate-derived mortality while questing for a host, whose densities depend on the rule of habitat perception. The relative densities of host types regulate the allocation of ticks to hosts. When host densities were higher in a given patch, the probability of tick-host encounters increased and the tick population had a lower mortality in the questing phase and an increase in the rate in which ticks progressed from stage to stage and reproduced.

Design concepts Emergence
Tick population dynamics are the results of the time of the year, the movements of the hosts through patches of suitable vegetation and the effect of vaccination. The time of the year has an impact because the natural phenology of the tick (as related to the weather) and the time after vaccination. Vaccination at a given period of the year will affect more or less feeding ticks on WTD. We tested the impact of an annual vaccination scheme conducted for 3 consecutive years in the complete study area on days 1, 70, 130, 190, 250 and 310. At each decadal, mortality was computed for ticks feeding on vaccinated WTD and dead females were removed from the tick population. Tick fertility was then calculated for surviving females according to the immunization time and the moment of the year, and both values were introduced into the model.

Adaptation
The adaptive traits in the model are the movements of WTD, which are driven by the carrying capacity of the patch that in turns depends upon its area and the distance to near patches in a spatial 2D network of movements.
Fitness. Reproduction of the ticks was performed according to simple rules that depend on the prevalent weather (see submodels below) and is affected by the vaccination protocols. For simplicity, we assumed that all hosts were adults in a stable population where changes in abundance, mortality and newborns did not occur. It was also assumed that all WTD were vaccinated in each simulation, to avoid the effects of multiple levels of immunity in a large population of moving WTD.

Interaction
There is no direct interaction between individuals, but indirect interaction through density dependent foraging.

Sensing
Movement of WTD relies only on physical attributes of the vegetation patch. We used the concept of habitat use by moving hosts, carrying and spreading ticks over a simulated territory (henceforth "landscape"). Two basic terms of this framework are traversability and recruitment [10]. Traversability is understood as permeability and addresses the importance of patch network features that drive the connectivity of the territory and therefore the movements of the animals. Recruitment focuses on the estimation of tick abundance at the patch level, linking the traversability with the tick's life cycle. The recruitment of a patch is a measure of its importance within the general network and is related to the abundance of ticks in that particular patch [3].

Stochasticity
The input data are deterministic.

Observation
The key output monitored from the model was the reduction of the cattle tick, R.

Initialization
The model of tick life cycle was run for 5 consecutive years at 10-days intervals.
After constant tick populations were obtained (e.g. each year has the same phonological component), we introduced parameters affecting tick reproductive rates by modifying the mortality rates and fertility of the feeding females on vaccinated hosts. The purpose was to simulate the effect of WTD vaccination because vaccination with tick protective antigens results in increasing tick mortality and reduction in the fertility of the surviving female ticks after feeding on vaccinated animals [11].
The host population was allocated to the vegetation patches, with a maximum carrying capacity of 5 WTD/ha and 10 cattle/ha [19,30,35,39]. At every decade (10days interval) of the year, hosts move over the network of patches according to their physical features (size and distance). Host movements are driven by equations governing the probability that a host move to another patch according to its size and the distance between the patches (see deer submodel below).

Input
The only external inputs were monthly temperatures and water vapor values (as saturation deficit) obtained from gridded climatology at 1 km resolution (available at http://www.climond.org). The original time resolution of the dataset (one month) was converted to 10 days interval by interpolation using splines.

Submodels Tick submodel
We used published results to produce a process-driven model calculating the phenology and abundance of the cattle ticks, R. annulatus and R. microplus, under field conditions. All the equations presented below have been developed, validated and explained before (see references for each entry). The tick life cycle was divided into several events or processes (e.g., oviposition, molting rates, mortality) that are driven by weather conditions. All the processes are included below, together with the complete explanation of the equations and the pertinent literature data. For every equation, T is the temperature in degrees Kelvin, and SD is the saturation deficit of the air in Hectopascals.
Rates of completion of oviposition per 10 days. Preoviposition decreases to a minimum of 2 days around 33ºC. The lower boundary is set at 12ºC [12].
Rates of female mortality per 10 days. A cohort of engorged females dies if it accumulates more than 83 mmHg in 10 days. To represent influence on female mortality of temperature alone, a degree-day model developed [13] was used.
The production of eggs by the female is a product of the temperature [14] according to the equation. This is the pure conversion of the female's body mass into eggs, and water contents of the air has not effects on it (but on the egg mortality rates, see below).
To estimate the time until first hatch of an egg cohort, a function of development rate, the rate-summation embryonic development model for R. annulatus [15] was used.
RHO25 represents the hourly development rate at 25ºC, TH the temperature at which rate-controlling enzyme becomes half active and half high-temperature inactive, TL the temperature at which rate-controlling enzyme becomes half active and half lowtemperature inactive, HH the change in heat content associated with high-temperature inactivation of the enzyme, HL the change in heat content associated with lowtemperature inactivation of the enzyme, HA the heat content of activation associated with the reaction catalyzed by a rate-controlling enzyme, r the universal gas constant, T the ground temperature.
Mortality of eggs depends upon the saturation deficit. Egg cohorts that complete development before the mortality proportion equals 1 begin to hatch, at which time the model applies the current mortality proportion to the cohort. [12,15].

(f) Duration of egg incubation
Published data [15] showed that mean duration of egg hatch at different temperatures equaled 77% of the time required to achieve first hatch. For surviving eggs in each cohort, the model approximates temporal distribution of egg hatch with a Weibull distribution [12].
where x is the normalized development time (multiples of time to first hatch) (g) Larval mortality rate [12,14].
Larval mortality rate is a function of temperature and water saturation deficit. Larval mortality rate is 1 in temperature < -5ºC or > 60ºC. Simulated larvae remain in the pasture until they die or attach to a host. Larval host-finding rates are explained in the next paragraph.
The rate at which questing larvae on vegetation infest cattle depends upon the complex and highly variable process of host-parasite contact. To estimate the hostfinding rate (HFR), which is the proportion of larvae picked up daily by a herd of cattle, we used equations that calculate the base host-finding rate (BHFR), a temperature-effect multiplier (TE), and a larval-density effect multiplier (LDE) [14].
The base hourly host-finding rate increases as host density (D, cows/ha or deer/ha) increases [12], according to the following formula: The multipliers used to modify this value, TE and LDE, have values bounded between 0 and 1. TE represents the increase in larval-questing activity prompted by an increase in temperature.
The on-host mortality rates of either R. microplus or R. annulatus are complex relationships depending upon the age of the hosts, their immune status and the tick loads. Further complications arise because we are modeling different mortality rates produced by two different hosts (deer and cattle). We used the mortality rates previously derived and validated [14,16]. We considered the same mortality rates for ticks derived from deer based on field results [11].

Deer submodel
White tailed deer is widely distributed in Mexico [1]. We primarily used the information regarding the expected abundance of WTD in the complete country as summarised in http://www.conabio.gob.mx/informacion/gis/?vns=gis_root/biodiv/distpot/dpmamif/d pmartio/odo_virggw (accessed on March 2012). This information was produced in the year 2009 based on published data [1]. Information about abundance and densities of WTD in the target area was updated with information obtained from 968 WTD ranches in the area of study. We assumed that WTD distribution follows the patches of adequate vegetation and that its density is proportional to the size of each We used Manifold 8.0 software (www.manifold.net) as a standard GIS package to download the spatial information, and produce the final landscape for further modeling. The capabilities of the GIS software produced the necessary data for calculation of distances, area of patches, and rules for movements of WTD. WTD move across the network of patches by rules governing movements according to the size of a patch of vegetation and its distance to another patches in the network. Thus, a patch must have a minimum area to be colonized; the larger the area, the more WTD will visit and stay; the shorter the distance to near patches, the less time WTD will remain at the same patch. Specific routines for creation of networks are available within the software package, which produced the complete picture of the landscape with routes of movements and spatial details of every patch. To calculate the details of the life cycle of the ticks (see below) a script in R [17] has been produced.
The host movement probability from patch i to patch j via frontier ij depends on the habitat perception rules. According to graph theory [10], the probability that an individual in node i will disperse to node j can be expressed in the form of a flux rate or dispersal probability matrix. Thus, the expected dispersal flux from patch i to j is: where Si is the area of patch i, S tot is the sum of the areas of every available patch, and p'ij is the probability of dispersal from i to j. This probability of dispersal p'ij is directly related to the area of patches i, j, and inversely related to the distance between them. Total traversability is thus defined as the sum of partial dispersal flux probabilities for every link, as a measure of the permeability of that patch to propagules coming from different patches in the network of the landscape [16,18]. To define p'ij we used a function, called habitat perception, of the form: where dij is the distance between the centroids of patches, r the mean dispersal distance of the host (km) and the gamma function. The parameter v relates the proportion of dispersing hosts as a response to the patch size. It has been assumed that proportion of moving hosts has a simple inverse relationship with patch size (i.e. small patches allow high migration rates) and the v parameter is simply a modifier of such a response in the gamma function [19].
The model for tick survival and growth rates runs in parallel with the host movement at 10-days intervals. Ticks find hosts depending on host density in the patch and the tick questing rates driven by climatic features. Therefore, the lower the density of WTD, the larger the time of quest and the higher the mortality driven by the weather and explained before in the equations for the tick submodel. Ticks have the same preference to feed on cattle than on deer [16] and the parasitic rates on hosts depend only on host availability at the patch. At each decadal, hosts move to another patch carrying a variable number of ticks according to the rules governing host-tick encounters. Engorged tick females drop from the hosts and colonize patches as "visited" by the hosts. After time, ticks are distributed over the network of patches visited by WTD and considering the effect of the climate, which describes tick development and mortality rates. Models were run separately for R. microplus and R.
annulatus because their performance is different under the same climatic features [14] and results were then summarized together.

Submodel of deer vaccination
The base mortality rates were further changed after the simulation of a vaccination event. Based on published data [11] we simulated the vaccination at a given moment, by increasing the mortality and reducing the fecundity of the females fed on vaccinated deer (no vaccination is simulated on cattle). The duration of immune response on the mortality has a temporal effect, and after peaking at a maximum it decreases slowly in time [11]. The effect is also different for either R. annulatus or R.
microplus. This is simulated by polynomial regressions, in which the effect of the time is evaluated in steps of 10 days. Figure 5 summarizes all the processes and submodels.
The equations for the increase of mortality for R. microplus (MORT_RM) and R.
annulatus (MORT_RA) are, respectively: The equations relating the decrease in fecundity with the vaccination are, respectively: