Time to Emergence of the Lyme Disease Pathogen in Habitats of the Northeastern U.S.A

TL;DR

This study estimates how long it takes for the Lyme disease pathogen to become established in new areas of New Hampshire, based on tick life cycles and host interactions.

Contribution

The study introduces a combined modeling approach to estimate the time to endemicity of Lyme disease in new habitats.

Findings

The time to reach Borrelia endemicity in New Hampshire ranges from 8 to 20 years.

Susceptible tick host populations strongly influence the time to endemicity.

Temperature, infection probabilities, and host populations are key factors identified through sensitivity analysis.

Abstract

Lyme disease is contracted by humans through the bite of a black-legged tick. The pathogen is maintained in tick populations through reservoirs of tick hosts in the wild, including pathogen hosts such as the white-footed mouse and mobile tick hosts such as the white-tailed deer. The pathogen may be introduced into a new location by the arrival of infected ticks, but it may take years to establish at endemic levels. A process-based model describes the tick life cycle, host interaction, and disease transmission. A statistical model predicts the populations of deer and mice across New Hampshire. We combine both types of models to estimate how long it takes this disease to establish in the environment and to identify major contributing factors.

Ticks carry a range of pathogens, the best known of which causes Lyme disease, prevalent in the northeastern United States. Emerging diseases do not yet consist of a wide range of Lyme diseases, raising the question of how long it takes for a newly introduced tick-borne disease to establish itself. The aim of this study was to address this question, with the agent of Lyme disease used as the test case. A prior process-based model of the Ixodes scapularis (Say 1821) life cycle and the transmission of Borrelia burgdorferi (Burgdorfer 1982) between this tick and its various hosts was used to predict the dynamics of disease introduction into a new area. The importance of temperature, infection probabilities, and tick host populations, relative to that of other factors, was established by a global sensitivity analysis using Latin hypercube sampling. The results of those samples were analyzed to determine the time to near-equilibrium. Eight locations in New Hampshire were chosen for high/low temperature, high/low mouse, and high/low deer values. Mammal abundance was estimated by relating the known mammal density from previous studies to a MaxEnt analysis output. The time required to reach Borrelia endemicity in the ticks of New Hampshire ranged from 8 to 20 years in regions where the tick population is viable, with a strong dependency on susceptible tick host populations.

Figures

Full text

1. Introduction

Lyme disease in humans was first identified in the northeastern U.S. in the 1970s due to a cluster of cases in Lyme, Connecticut. Its agent, the bacteria Borrelia burgdorferi (Burgdorfer 1982) [1], and its vector, the tick Ixodes scapularis (Say 1821) [2], were identified in 1981 [3]. From 1985 to 2013, reported cases spread northward, with few initially reported in New Hampshire [4,5,6]. Lyme disease is now the most common arthropod-borne disease in the U.S. [3]. New Hampshire has one of the highest rates of Lyme disease in the United States, with it having been identified in all 10 counties, according to the New Hampshire Department of Health and Human Services. The ticks themselves are present throughout New Hampshire but are more prevalent in southern areas, according to recent passive surveillance [7].

The means of transmission is between ticks and susceptible reservoir hosts [8]. The geographic expansion of disease risk is likely due to tick hosts that are mobile, such as deer and ground-nesting birds [8,9]. Deer do not host the bacterial agent but may carry ticks that are already infected [8].

The vector I. scapularis, or the black-legged tick, is a hard-bodied tick that requires up to three blood meals to complete its life cycle. Four stages (egg, larva, nymph, and adult) are punctuated by blood meals obtained by questing behaviors that expose the tick to weather. The ability of the tick to mature to the next stage depends on its success in finding a host, which, in turn, depends on the abundance of and access to hosts [10,11]. Intermediate periods involve digestion and molting, which take place underground or in leaf litter [8,12]. The length of these intermediate periods depends on the temperature. Humidity is also a factor [13]. Therefore, the progression of Lyme disease across the landscape depends on the underlying ecosystem supporting the tick population.

A process-based model was used to simulate tick and host populations and disease transmission over time. A sensitivity analysis was conducted to confirm the importance of temperature and host distribution to tick populations and infection prevalence in ticks. A data-driven statistical model was used to estimate the abundance of two important tick hosts in New Hampshire. Eight locations were selected for further study and comparison using the process-based model.

This study shows the importance of temperature in determining the viability of tick populations and the importance of host distribution in determining the time to steady state of the percent of nymphal ticks infected, at which point the disease is completely endemic in the tick population. In addition, it is shown that the time to the endemicity of disease in the tick population ranges from 5 to 20 years, with a high dependence on the population of mice, which are competent Borrelia hosts, and a low dependence on the population of deer, which are not Borrelia-competent hosts but play a large role in feeding and transporting ticks to new locations.

2. Materials and Methods

This study adapted a prior model, which is described here, with the equations given in the Appendix A. This model was based on the life cycle, temperature dependencies, and host requirements of I. scapularis, tuned to the temperature profile of Hanover, NH, and adjusted with the day length-related diapause to match the qualitative characteristics observed in previous studies [13,14,15]. The model was the basis of a global sensitivity analysis, which established the importance of several parameters, including some involving tick hosts.

Two of the main tick hosts are mice and deer. A statistical model of the relative abundance of these hosts was developed for the state of New Hampshire. Many researchers have examined habitat use/selection in accordance with niche quantification by measuring the environmental parameters of habitat types [16,17,18,19,20,21,22]. The determination of factors and the relationship of multiple variables is appropriate for habitat use studies [23]. The studied habitat attributes of an organism are associated with the “realized” niche, since observational field data merely examine a possible niche volume and not the entire range that a species may exist in. An understanding of the habitat can greatly aid conservation efforts and ecological studies [24,25].

Some robust analytical methods use presence-only data, sampling from background data, to develop a statistical model that extrapolates across the landscape habitat probability [26]. One such method is maximum entropy modeling. Using MaxEnt software version 3.4.1 [27], we estimated a site probability map across the landscape. Researchers have used this habitat estimation method on regional to landscape scales to infer archaeological settlement patterns in the Brazilian Amazon (black earths) [28], Bolivian bamboo forests (geoglyphs) [29], and Michigan [30]. This modeling approach has also been used in Southern Maine to estimate the habitat probability of New England Cottontails [31].

Based on this analysis, eight locations were chosen for further analysis that vary by the temperature profile, habitat probability of susceptible rodent hosts, and habitat probability of disease-resistant deer hosts.

2.1. Model Overview

The process-based model used in this study was based on prior models that include the life cycle and disease dynamics of both the arthropod vector I. scapularis and its hosts [14,15,32]. This tick has a multi-stage life cycle from egg to adult, with three molts, three questing periods, and three blood meals. Temperature-dependent maturation rates control parts of the tick life cycle occurring off-host, while the availability of hosts determines on-host survival rates [13]. In addition, ticks have inactive periods in which they do not quest for a host; these appear to be unrelated to temperature and may be related to day length [33]. These periods are sometimes referred to as “diapause” and have also been incorporated into prior models to explain observed questing patterns in the field.

Six classes of host were included, categorized as competent (able to carry and transmit B. burgdorferi) or incompetent and mobile or stationary. Competent hosts were further sorted into infectious or susceptible hosts. Competent stationary hosts are small rodents, in particular, the white-footed mouse (Peromyscus leucopus (Thomas 1895) [34]), which is the primary host of larval ticks. Competent mobile hosts are ground-nesting birds. Incompetent stationary hosts include mid-size animals, such as opossum. Incompetent mobile hosts are exclusively white-tailed deer (Odocoileus virginianus (Zimmerman 1780) [35]), which, while unable to transmit disease, can transport an infected nymph for a considerable distance. Each of these hosts has simple population dynamics limited by an environmental carrying capacity, as well as a per-individual capacity for hosting ticks [11]. Equations and default model parameters are given in Appendix A.1. All parameters are based on prior work with this model, in which multiple sources of data in the northeast were used either to estimate a parameter directly or to adjust parameters to produce outputs comparable to those observed in field studies [14].

2.2. Sensitivity Analysis

The model for the tick population and disease dynamics has 71 parameters describing temperature effects, host dynamics, diapause due to day length, the efficiency of questing for hosts, and the carrying capacity for both hosts and ticks on hosts. A global sensitivity analysis was performed using MatLab’s Latin hypercube sampling code, with all 71 parameters, repeated 100 times [36]. The specified range for each parameter was ±20% of the default parameters, as shown in Table A1 in Appendix A.2. The initial conditions were drawn from prior work.

The model describes 1 km^2^ in which ticks are already present but not yet infected. The initial conditions were set for 1 January of year 1, with 6,453,100 eggs, 4,856,100 engorged uninfected larvae, and 291,360 engorged uninfected nymphs. Disease was introduced by assuming that 100 infected nymphs were present on a deer. The remaining tick categories were all set to zero. The number of eggs, percent of nymphs infected, and percent of adults infected on 1 January of year 46 were calculated for each of the 100 runs, and the correlation coefficient was obtained for each of the 71 variables calculated.

Although there was a large annual variation in populations, the value of these populations on 1 January reached equilibrium, as did the infection prevalence in nymphs and adults. The approximate time to 1 January equilibrium was calculated as the first year in which the values were within a given threshold of the final year’s value, for each of the 100 runs specified in a hypercube sample. As the number of eggs was large, the threshold for eggs was set to 50,000. The other two quantities are given as percentages, with a threshold of 0.02 or within 2% of the final value for the “percent infected” nymphs and adults. All measures were taken on Jan 1 of each year, as only a few stages are present at that point. Time series were graphed for the three quantities of interest for all 100 runs to visually confirm calculations.

2.3. MaxEnt for Host Distributions

Statistical models were used to determine the locations of high and low deer or mouse densities. Spatial modeling was performed using Maximum Entropy Species Distribution Modeling Software 3.3 [26]. Geospatially derived layers were used to calculate the driving factors for site location and to estimate the probability of sites being suitable for a given species across the research domain, in essence a habitat probability map. The use of maximum entropy modeling allowed us to characterize the driving environmental and spatial indices for multiple species. Species point information was used for two species of host vectors (white-footed mouse and white-tailed deer).

The source locations for deer and mouse were obtained from the Global Biodiversity Information Facility (GBIF), an international network and data infrastructure funded by the world’s governments and aimed at providing anyone, anywhere with open access to data about all types of life on Earth. For mice, 1704 location/occurrences in NH were used. One study was removed from the analysis that was very localized. For deer, 1418 location/occurrences in NH were used.

The MaxEnt model was run with cross-validations when appropriate. Iterations were run using a convergence threshold set to 0.0001, with a default prevalence of 0.5. Response curves were generated for each input variable in our MaxEnt model. Jackknife tests were run to examine the importance, or relative contribution, of each variable to the model. The jackknife tests allowed for an examination of how a given variable contributed to the overall model, both when it was used exclusively to build the model and when it was excluded from the model. The jackknife tests were performed for both training and test datasets. The area under the curve (AUC) statistic was used to gauge the predictive capacity of the model and how well it performed compared with a null model. Models with AUC values that are greater than 0.75 are considered to predict the test point distribution accurately [30].

Geospatial layers provide insights into the distribution and driving attributes of ecological phenomena, such as species distribution. A wide range of variables that were likely to be insightful and exploratory were utilized. The following layers are described below, and their inclusion in specific model ensemble runs are mentioned in each of the specific sections. Table A2 describes each dataset and its spatial resolution. Table A3 and Table A4 give estimates of the relative contributions of the environmental variables to the MaxEnt model for P. leucopus and O. virginianus, respectively. The MaxEnt model was run at a 1 m resolution across the entire state of New Hampshire. Imagery or geospatial datasets that had a lower resolution were resampled to 1 m for modeling purposes.

2.4. Specific Locations

Based on the results of the MaxEnt model and these measurements, eight locations were chosen in New Hampshire with either high or low temperature profiles, high or low P. leucopus densities, and high or low O. virginianus densities. The MaxEnt model was used to find locations with high and low precipitation and temperature, as well as high and low probabilities of mouse and deer habitats. These locations were then used to extract the time series of meteorological time series data to drive the tick model. Temperature data were fit with a Fourier series to parameterize the model for each location. The temperature data were obtained from Daymet [37]. The 8 locations were not field-based or sampling locations. These locations were used to extract long-term meteorological data on temperature and precipitation to drive the model.

Recent population density estimates from around the state range from generally less than 4 deer per square mile of habitat in the White Mountains to 18–24 per square mile in more southern locations, and, as mentioned above, they have been averaging at about 12 per square mile statewide [38]. For low deer populations, 1.54 deer per km^2^ (4 per square mile) was used, and for high deer populations, 7.72 per km^2^ (20 per square mile) was used.

Although multiple rodents and other small mammals may carry B. burgdorferi, the white-footed mouse is an important reservoir species for P. leucopus [8]. Donnelly et al. found population densities of 5 to 39 per hectare, or 500 to 3900 per km^2^ [39]. A separate study by Wolff et al. found that home ranges averaged at 590 m^2^, with a minimum size of 500 m^2^, consistent with Donnelly et al. [40].

3. Results

3.1. Model Performance

Tick populations have large seasonal effects, with questing adults present early and late in the season and questing nymphs found primarily in the middle of summer. As the Lyme disease pathogen advances in tick and host populations, little seasonal effect was seen in the model on the percent infected nymphs or adults. The final percent infected nymphs ranged widely among the 100 runs, from 6 to 20%, while the final percent infected adults ranged from 10 to 20%.

3.2. Sensitivity Analysis

Egg populations were measured on 1 January of each year of the model to suppress seasonal effects. Equilibrium was reached in every run, which can be visually checked in Figure 1A. The final number of eggs was 0.5 × 10^6^ to 2 × 10^6^ across all runs, as shown in Figure 1D. The parameters with the strongest positive correlation with egg abundance were the birth rate (b), the transition of feeding larvae to engorged nymphs (m3c), and the probability that a questing adult finds any host (qA). Those with the strongest negative correlation with egg abundance were the death rate of eggs (de), death rate of engorged nymphs (dn1), and death rate of engorged adults (dA4). However, no single correlation was large, and many parameters had a similar effect size, as seen in Figure 1B. The time to egg equilibrium was fairly short, as seen in Figure 1C.

The percent infected nymphs was measured on 1 January and can be seen to reach equilibrium in Figure 2A. The parameters with the strongest positive correlation with the percent of nymphs infected at the end of the run were the mean annual temperature (tempMean), the probability of a feeding larva picking up infection from an infected host (pL), the probability of a competent uninfected stationary host (i.e., a rodent) picking up infection from an infected feeding tick (pCUS), the temperature dependence of the nymph maturation rate (tmfn), the death rate of competent uninfected stationary hosts (i.e., rodents) (dCIS), and the carrying capacity of competent mobile hosts (i.e., birds) (KCM). In this case, the parameters with the strongest correlation stood out a bit more than those for egg populations, as seen in Figure 2B. The range of times to equilibrium and the infection prevalence in nymphs at equilibrium are shown in Figure 2C and Figure 2D, respectively.

The percent infected adults was similar, with a year longer time frame due to the extra year of maturation, as seen in Figure 3A,C. The parameters with the strongest positive or negative correlation with the percent infected adults were the same as those for nymphs, as seen in Figure 3B. The final infection rates were a bit higher than those for nymphs, as seen in Figure 3D.

3.3. Years to Equilibrium and Endemicity of Disease

The initial conditions were the same for all runs, but the equilibrium values for the final egg populations varied with the parameter shifts in the 100 runs shown in Figure 1D. Figure 1C shows a histogram of the time to equilibrium, with the majority of runs approaching near-equilibrium in 4 years or less. In contrast, the number of years to equilibrium for the percent infected nymphs ranged from 8 to 21 years, as seen in Figure 2C, while the number of years to equilibrium for the percent infected adults was, predictably, a year longer, ranging from 10 to 22 years, as seen in Figure 3C. The sensitivity analysis suggests that temperature and host abundance strongly affect the time to steady state.

3.4. MaxEnt for Host Distributions

Figure 4 and Figure 5 show the results of the analysis of the mouse and deer distributions in New Hampshire. Figure 4a and Figure 5a show the sensitivity versus specificity for the mouse and deer predictions. In both cases, the AUC was over 0.75, indicating a good prediction. The MaxEnt model provides an indication of the key habitat of the modeled species, and this is often associated with a higher population density. Note that areas of high/low habitat probability are not particularly correlated between mouse and deer. This allows for a selection of sites with a high/low temperature, mouse habitat probability, and deer habitat probability.

3.5. Comparison of Eight Locations

The eight locations in New Hampshire chosen for comparison varied by temperature and host density. Densities were set to high and low values found in the literature. The four locations with low temperature profiles could not sustain the tick population. The four locations with warm temperature profiles were able to sustain tick populations. These are summarized in Figure 6. The time series for one of these is shown in Figure 6B, with seasonality. The sites with higher mouse populations reached a steady state earlier than those with low mouse populations, as shown in Figure 6A. Deer populations did not seem to make a large difference.

4. Discussion

This study makes use of a model previously described in the literature [14,15,32]. The maturation rates of I. scapularis in the model are temperature-dependent and parameterized based on experiments by Ogden et al. [13]. These dependencies prevent efficient maturation in regions with short warm seasons. The strong link between disease prevalence and temperature is confirmed in this study by the sensitivity analysis, which indicates that the mean annual temperature is one of the parameters to which the B. burgdorferi prevalence in nymphs was most sensitive, consistent with other studies [14,32,41,42,43,44]. In the four specific locations modeled that had a low mean annual temperature, the model showed tick populations declining to zero over time. Two caveats must be mentioned about this result. First, it may well be that some microclimates in a region may be quite cold, while others nearby could be warm enough to sustain tick populations. Therefore, it would not necessarily be true that there are no ticks in a given town or square kilometer of land. Second, ticks spend much of their maturation time hiding in soil or leaf litter, or even in animal dens, which have warmer temperatures than the air during the winter [8]. Measurements of these subsoil habitats would be useful in improving the temperature dependence of the model. It is worth noting, however, that passive tick surveys do not report ticks from colder regions of New Hampshire [7].

The strong links between Lyme disease and host populations are also confirmed by the sensitivity analysis of the model. Two of the most influential parameters were the death rate of competent uninfected stationary hosts (i.e., rodents, such as mice) and the environmental carrying capacity of Borrelia-competent mobile hosts (i.e., ground-nesting birds). The default parameters for the sensitivity analysis included the populations of these two categories, as well as two types of Borrelia-incompetent hosts. The difficulty of estimating animal populations across a region was addressed through MaxEnt software, which was able to estimate the relative abundance of the Borrelia-competent white-footed mouse and the Borrelia-incompetent deer populations across the state of New Hampshire. For the four warm-temperature sites, the rate of disease establishment at a given location had a strong dependency on mouse populations but no observable dependency on deer populations. This dependency only concerns the time to pathogen equilibrium and the prevalence of the pathogen in ticks. It does not reflect the number of ticks present, which certainly depends upon all hosts, nor does it represent the number of infected ticks, which determines the disease risk for humans.

Based on the sensitivity analysis, the range of time to disease endemicity in nymphs varied from 5 to 20 years, depending on the parameters, with 12 years being the mode of the distribution. This finding is consistent with the reported spread of Lyme disease in New Hampshire between 1998 and 2019, which increased dramatically in 20 years [5]. Because deers are highly mobile transporters of ticks, infected or not, it is likely that disease enters a town from a neighboring town. The time lag between disease reported in a county adjacent to one that has already had cases is estimated to be 7 years [45].

The analysis of the four warm-temperature sites showed the model to be at near-equilibrium in 10 years for high mouse populations and at equilibrium at 15 years for low mouse populations, consistent with the sensitivity analysis. As Borrelia-competent hosts, mice are believed to play a major role in the transmission of the Lyme disease spirochete [8,11,46,47]. Both the sensitivity analysis and the model results for warmer sites confirm this conclusion.

Sites with warmer temperatures were chosen in the towns of Plymouth (low mice and low deer populations), Madison (high mice and low deer populations), Deerfield (low mice and high deer populations), and Jackson (high mice and high deer populations). Jackson had 78.3 and 59.2 cases of Lyme disease per 100,000 in two five-year intervals (2005–2009 and 2010–2014). Madison had rates of 33.8 and 117.7 in the same intervals. Plymouth (14.3 and 62.4) and Deerfield (74.7 and 225.3) were the sites with lower mouse populations [48]. It is clear from these data that the eventual prevalence of disease in the human population is not particularly related to the time frame of disease establishment in the tick population. In other words, this study does not measure the risk of disease for humans. That risk is present long before the tick infection rate reaches a steady state, and it is dependent on many factors. These include human behaviors such as picking up the pathogen in one area and reporting the associated disease in another, the time spent outdoors, and other behaviors. Additionally, areas with low mice populations may offer an abundance of other disease reservoirs not included in this model.

5. Conclusions

In this study, we addressed the question of how long it takes a newly introduced disease to establish in a new location, using a process-based model of the tick life cycle and data-driven estimates of host populations. We found a time range from 5 to 20 years, depending on the tick and disease host populations and temperature. Model tick populations were not sustained at very low temperature ranges. In warmer temperature ranges, the disease prevalence in ticks was highly dependent on mouse (or other competent host) populations.

Although based on the transmission of B. burgdorferi, this model also has implications for other tick-borne pathogens, such as those responsible for anaplasmosis (Anaplasma phagocytophilum) and babesiosis (Babesia microti) [49]. Although currently less prevalent, one may expect these pathogens to also spread geographically. Based on the case of B. burgdorferi, this study may provide insights into the future spread of other pathogens by I. scapularis.

Most aspects of the model involve the life cycle, hosts, and temperature dependencies of I. scapularis, which is the same vector for other pathogens. The main difference in modeling other diseases transmitted by I. scapularis is the variation in the transmission rates between the tick and the host and the demographics of competent versus incompetent hosts for a different disease. The model for B. burgdorferi transmission is highly sensitive to the transmission rates between the vector and the host. This result indicates that measuring vector–host transmission rates for emerging tick-borne diseases should be a high priority.