Journal of The Royal Society Interface
You have accessResearch article

Modelling the seasonality of Lyme disease risk and the potential impacts of a warming climate within the heterogeneous landscapes of Scotland

Sen Li

Sen Li

Environmental Change Institute, University of Oxford, South Parks Road, Oxford OX1 3QY, UK

[email protected]

Google Scholar

Find this author on PubMed

,
Lucy Gilbert

Lucy Gilbert

The James Hutton Institute, Craigiebuckler, Aberdeen AB15 8QH, UK

Google Scholar

Find this author on PubMed

,
Paula A. Harrison

Paula A. Harrison

Environmental Change Institute, University of Oxford, South Parks Road, Oxford OX1 3QY, UK

Centre for Ecology and Hydrology, Library Avenue, Lancaster LA1 4AP, UK

Google Scholar

Find this author on PubMed

and
Mark D. A. Rounsevell

Mark D. A. Rounsevell

School of GeoSciences, University of Edinburgh, Drummond Street, Edinburgh EH8 9XP, UK

Google Scholar

Find this author on PubMed

    Abstract

    Lyme disease is the most prevalent vector-borne disease in the temperate Northern Hemisphere. The abundance of infected nymphal ticks is commonly used as a Lyme disease risk indicator. Temperature can influence the dynamics of disease by shaping the activity and development of ticks and, hence, altering the contact pattern and pathogen transmission between ticks and their host animals. A mechanistic, agent-based model was developed to study the temperature-driven seasonality of Ixodes ricinus ticks and transmission of Borrelia burgdorferi sensu lato across mainland Scotland. Based on 12-year averaged temperature surfaces, our model predicted that Lyme disease risk currently peaks in autumn, approximately six weeks after the temperature peak. The risk was predicted to decrease with increasing altitude. Increases in temperature were predicted to prolong the duration of the tick questing season and expand the risk area to higher altitudinal and latitudinal regions. These predicted impacts on tick population ecology may be expected to lead to greater tick–host contacts under climate warming and, hence, greater risks of pathogen transmission. The model is useful in improving understanding of the spatial determinants and system mechanisms of Lyme disease pathogen transmission and its sensitivity to temperature changes.

    1. Introduction

    Lyme disease (or Lyme borreliosis) is the most prevalent vector-borne disease in temperate regions of the Northern Hemisphere [1]. In the UK, the annual number of confirmed Lyme disease cases is over 1000 and still increasing in some areas each year [2]. The causative agents of Lyme disease belong to a species complex of spirochaete bacteria named Borrelia burgdorferi sensu lato (s.l.), which are transmitted principally by the tick Ixodes ricinus (castor bean, deer or sheep tick). In Scotland, B. afzelii (associated with rodents) is the most abundant genospecies, followed by B. garinii (associated with birds) and B. burgdorferi sensu stricto (s.s.) (associated with both rodents and birds) [3]. As no vaccine is available, the prevention of Lyme disease relies heavily on improving understanding of the mechanisms of disease transmission and integrating multi-disciplinary knowledge and data to predict the patterns of disease risk in a changing environment.

    The risk of Lyme disease is strongly related to the density of active infected ticks, which can be influenced by a wide range of biophysical factors related to host communities, climate and landscape (see recent reviews [49]). Temperature, in particular, has gained much attention in eco-epidemiological studies of Lyme disease, in line with increasing concern about the effects of climate warming over the last few decades. Temperature affects the abundance of infected ticks in two main ways. First, it directly affects the behaviour (e.g. activity and diapause), interstadial development rate, fertility and survival of ticks [10,11]. Second, it influences the population and habitat suitability of host species. In addition, a warmer climate at higher altitudes may result in an increase in the population of important tick hosts such as roe deer, Capreolus capreolus [12]. Population dynamics of rodents (the transmission host) may also be influenced by temperature, as climate change shifts local faunal abundance and diversity [13].

    Mechanistic models, that explicitly represent the key population groups and processes underpinning their interactions, are capable of providing unique insights in support of health-related decision-making [14]. They have increasingly been applied to explain disease transmission system dynamics, understand observed disease patterns, test scenarios for potential future risk patterns and propose data and methodological improvements [15]. Existing mechanistic models for studying the effects of temperature on Lyme disease risk are increasing in number (e.g. [1620]). However, existing models do not incorporate spatial heterogeneity or focus on issues related only to the dynamics of vector populations, ignoring key components, such as host/pathogen distribution and habitat suitability, which in reality are distributed unevenly and can change rapidly over time. To bridge these gaps, we developed a mechanistic model to investigate the spatio-temporal dynamics of tick-borne diseases by integrating a number of temperature-dependent functions for the temporal dynamics of tick populations into a recently developed spatial model for Lyme disease risk. The model considers the heterogeneity in tick population ecology, host movements and pathogen transmission across the landscape [21,22]. The model was developed for Lyme disease risk assessment in mainland Scotland where the potential for disease exposure is high due to substantial public participation in nature-based activities.

    The principal objective of this study was to use the outputs of the mechanistic model for the spatio-temporal pattern of infected tick abundance as a biophysical risk indicator (or ‘hazard’) of Lyme disease. The model was designed to take account of current scientific understanding of the system mechanisms so that it can serve as an explanatory tool for further research and educational purposes. A secondary objective was to predict the disease risk pattern and its seasonality for mainland Scotland. We did not attempt to make precise annual predictions, nor predictions of the distributional spread of the pathogen, but rather focus on the long-term (i.e. decadal) average approximation of the seasonality, which is appropriate for application to scenarios of long-term environmental change. To gain greater insight into the spatial determinants of disease risk pattern and its sensitivity to temperature changes, two modelling experiments were conducted: (i) an exploration of the relationships between the predicted disease risk patterns and environmental factors of local interest, i.e. annual average temperature (affects tick survival, development and activity), elevation (modifies the temperature, habitat distribution and deer movement) and deer density (influences tick reproduction); and (ii) prediction of changes in disease risk pattern under simple climate warming scenarios (i.e. regional increases of temperature by 1°C, 2°C and 3°C throughout the year). Finally, we assessed the data and methodological barriers of adopting such a complex zoonotic disease model for regional-level real-world applications.

    2. Material and methods

    A mechanistic, agent-based model was developed by integrating recent developments in simulating the temperature-driven temporal [1620] and spatial [21,22] dynamics of I. ricinus populations and B. burgdorferi s.l. transmission in mainland Scotland. Multi-sectoral data were used to prepare model inputs [2330], modify model functions [23] and evaluate the model outputs [24] (see details and discussions in the electronic supplementary material). The model was used to predict the spatio-temporal pattern for the density of infected nymphs (DIN, as a Lyme disease risk indicator) based on 12-year (2000–2011) mean weekly temperature surfaces. The model was then applied to a number of simple temperature warming scenarios to explore the possible consequences of climate change on the DIN pattern. Finally, we further examined the spatio-temporal patterns of simulation predictions by summarizing the predicted DIN changes according to annual average temperature, elevation and deer density parameters to assess how the effects of these parameters may vary with a warming climate.

    2.1. Model overview

    The conceptual framework for the Lyme disease risk model is presented in figure 1. It is programmed using Repast Simphony (v. 2.2) [31] in which the environment is represented as a two-dimensional, rectilinear grid with a cell size of 1 km2 and a time step of one week. There are three interactive layers within the model where cells represent the population distributions of ticks and host animals and the configuration of the landscape. The mechanisms for the spatio-temporal developments and interactions of the model layers are represented as transition rules, which are grouped into three sets concerning tick population dynamics, pathogen transmission, and host population dynamics and movement patterns, respectively. The framework was designed to describe the general ecological processes of Lyme disease but the present model was parametrized for mainland Scotland. A summary of the model is provided in this section and detailed descriptions of the layers, parameters and transition rules are presented and discussed in the electronic supplementary material (§S2, see all parameter values in table S2.1).

    Figure 1.

    Figure 1. Model conceptual framework: interactions between ticks, hosts, pathogens and landscape. (Online version in colour.)

    2.1.1. Tick population layer

    Four life stages are modelled: ‘egg’, ‘larva’, ‘nymph’ and ‘adult’. In each post-egg life stage, ticks could be in questing, feeding or interstadial development phases. Total and infectious populations in all stages are simulated for each cell. When encountering hosts, questing ticks attach for blood meals, then drop and develop into the next life stages. Female adult ticks (assuming half of the emerged adult ticks are females) also need blood meals to produce eggs that hatch into larval ticks. Durations of feeding are assumed to be less than one week for larvae, and one week for nymphs and adults. We also assume that questing activity and interstadial development are sensitive to temperature, while feeding success is assumed to be dependent on the density of hosts [1620].

    2.1.2. Host population layer

    Populations of three generalized host types are simulated for each cell: (i) transmission hosts that are capable of transmitting B. burgdorferi s.l. (e.g. rodents [32], birds [33] and lagomorphs [34]), (ii) deer (e.g. roe deer and red deer) and (iii) livestock (e.g. sheep and cattle). In the model, the overall number of transmission hosts varies with season, but that of deer and livestock is fixed. The spatial distributions of hosts can vary between time steps due to animal weekly and seasonal movements (i.e. home-ranging and dispersal), resulting in ticks being transported from one place to another. Owing to limited data, only the generic B. burgdorferi s.l. was considered in the model, not the individual genospecies. As the density of common rodent species [35] is much greater than that of birds [30] in Scotland, the transmission host demography is assumed to be largely dependent on rodents. Hence, the populations of bird and lagomorph are included in the host layer, but their dynamics are not modelled explicitly. Therefore, considering the scale (the whole of mainland Scotland) and spatial resolution (1 km2), this model is more relevant for rodent-specialized Borrelia genospecies such as B. afzelii and B. burgdorferi s.s. than for the bird-specialized B. garinii. Moreover, the movement of transmission hosts is assumed negligible for a cell size of 1 km2, as the home range sizes of rodents and birds are usually around 0.1–0.2 ha, and the dispersal distances are less than 1 km for rodents and at the continental or global level for birds [36,37]. In the model, deer and livestock are assumed to be tick reproduction hosts which provide a large quantity of blood meals for tick development and reproduction but which do not support systemic B. burgdorferi s.l. transmission [3841]. Sheep may support non-systemic pathogen transmission between co-feeding ticks [42]. However, such a transmission route is not included in the present model as there is debate on its significance in the maintenance of B. burgdorferi s.l. [43,44]. The two reproduction host types (deer and livestock) use different land covers and perform different movement patterns. Home range sizes of both deer and sheep are assumed to change seasonally. The seasonal dispersal of deer (i.e. downhill migration during winter and uphill migration during summer [45]) is also modelled explicitly. While ticks in each post-egg life stage are assumed to be capable of feeding on all host types, in the model larvae preferred transmission hosts (rodents and birds), adults preferred reproduction hosts (deer and sheep) and nymphs were assumed to be generalists.

    2.1.3. Landscape layer

    The landscape is represented by four cell states: woodland, grassland, heathland and non-vegetated/other areas. The model assumes that ticks can inhabit grassland and heathland, but that woodland is more suitable [5,12]. Mortality rates of ticks, densities of hosts and movement patterns of hosts vary with the four land cover types. Transmission hosts can inhabit all woodland, heathland and grassland. Deer are considered to inhabit mainly woodland and heathland, and they can spend a proportion of their time in grassland for grazing. Conversely, livestock are assumed to be present mainly in grassland, but are also able to graze in heathland and woodland. All hosts can enter non-vegetated areas, but do not stay, hence, it is assumed that no ticks drop off in such areas. Elevation is used to simulate the seasonal uphill/downhill dispersal of deer and altitudinal changes in their home range.

    At each successive time step, these transition rules are applied in sequence to update cell states simultaneously. Firstly, within each cell, the infectious and total populations of questing ticks at each life stage are updated by adding in the unfed ticks moulted from a previous life stage and removing those that die in the previous time step. At the same time, the infectious and total populations of transmission hosts are updated according to the seasonal carrying capacity and birth/mortality rates. Secondly, hosts encounter questing ticks in the same cell and pathogen transmission takes place. Finally, between-cell host movements are considered, i.e. all nymphs and adults attached to out-moving animals are transported and dropped off at the end of the time step. The transport of larvae, however, is considered negligible, as their feeding duration is assumed to be less than one week (around 3 days) and they are assumed to feed preferentially on small transmission hosts with small home ranges which do not move between cells.

    2.2. Model inputs, outputs and initialization

    2.2.1. Input

    Multi-sectoral datasets were processed and integrated, including empirical data on tick behaviour, B. burgdorferi s.l. infection and distribution [23,24], earth observations on land cover and habitat types [25], long-term publically accessible meteorological data [26] and host habitat suitability and distribution estimations based on census data, literature-based qualitative information and model predictions [2730]. Data sources and preparation details are presented in the electronic supplementary material (§S3).

    2.2.2. Output

    The model can output an infected/uninfected tick density surface for all the different life stages (egg, larva, nymph and adult) at difference phases (questing, feeding or interstadial development). In this study, only the density of questing infected nymphs (DIN) was output, as this is a widely used Lyme disease risk indictor [49]. Nymphs are an order of magnitude more common than adult ticks. Nymphs are also smaller than adult ticks and less likely to be spotted, allowing them to complete their blood meal and transmit Borrelia pathogens. Nymphs have been shown to be the stage causing most Lyme disease cases in people [46,47]. Furthermore, we focused on simulations in woodlands, as field studies have shown that many more ticks have been observed in forests than open habitats [48].

    2.2.3. Initialization

    All simulations were initialized with assumed initial densities of 2 × 105 and 5 × 104 per km2 for the total and infected questing nymphs, respectively, in woodland. All results were recorded after 2600 time steps (50 years) to ensure stabilized yearly cycles had been reached.

    2.3. Model evaluation

    Two types of model evaluation were undertaken. The first qualitatively compared the shape of the overall tick/pathogen seasonality curve with empirical findings in Scotland or the UK from the literature. The second directly compared the spatial pattern of DIN with field observations based on: (i) species distribution data from the NBN (National Biodiversity Network) Gateway (http://data.nbn.org.uk/) to obtain a general sense of the degree of overall agreement and detect any under- or non-detected patterns and (ii) the field study on the infection prevalence in nymphs (NIP) [24] to check the model's predictive power at 24 forest sites in mainland Scotland.

    2.4. Model simulation

    The main objective of the simulation exercise was to predict the possible long-term disease risk pattern and its seasonality for mainland Scotland, as well as possible changes under climate warming scenarios. To gain an overview of disease risk pattern and seasonality, weekly surfaces of infected nymph density (a proxy of Lyme disease risk) were produced for mainland Scotland using a set of 12-year (2000–2011) average weekly temperature surfaces. Predicted disease risk pattern and its seasonality were compared between administrative regions. Then, the potential effects of climate warming scenarios were tested by applying increases of 1°C, 2°C and 3°C to the 12-year average temperature surfaces. Finally, by analysing the predicted DIN, we investigated the effects of some selected cell-level variables on disease risk, including: (i) annual average temperature (influences tick survival, development and activity); (ii) elevation (modifies the temperature, habitat distribution and deer movement) and (iii) deer density (affects tick reproduction). Temperature was selected due to climate change concerns [2,7,49], while deer densities were selected as these can be managed as part of disease risk controlling strategies [22,50]. Elevation was hypothesized to be able to capture the intricacies of temperature, hosts and habitat better than the other two variables [51]. We first checked that the three variables were only weakly associated with one another (Spearman's correlation coefficients were −0.31 between (i) and (ii), 0.39 between (ii) and (iii), and −0.18 between (i) and (iii), and all significant at the 0.05 level), before summarizing and plotting the predicted disease risk patterns against the three selected cell-level variables to determine the efficacy in predicting the effects of climate warming. To improve the visualization of plots, we reduced the size of the predicted DIN data by averaging the cell-level predictions for each interval of 1(m) for elevation, 0.1(°C) for temperature and 0.1 (heads/km2) for deer density.

    3. Results

    3.1. Model performance

    3.1.1. Stabilized yearly cycles

    Very similar stabilized yearly patterns in the DIN (less than 0.01% overall difference) were found after varying initial values by ±90%. Hence, the model outcomes seem to be largely independent of differences in the initial tick densities. Model predictions on the seasonal patterns of questing tick populations and infection prevalence in different habitat types are provided and compared with field evidence found in the literature in the electronic supplementary material (§S4.1). Stabilized yearly cycles of infectious tick populations were achieved after 20–40 simulated years in all simulations in woodland habitats. Tick populations in other habitats also achieved stabilized cycles at lower levels, even though no ticks were assumed to be present at initialization. The questing larval and nymphal ticks had similar predicted seasonal patterns with a peak in autumn (September–October). The questing adult females had a relatively symmetrical predicted pattern with peak values in summer (June–July). The nymphal infection prevalence (NIP) and adult infection prevalence were predicted to be relatively stable over the course of the year. To assess the overall model performance, additional results on the predicted transmission host populations and feeding tick populations are also provided and discussed in the electronic supplementary material.

    3.1.2. Spatial dynamics

    Predictions and discussions on the spatial patterns of tick and NIP are provided in the electronic supplementary material (§S4.2). In general, model performance in predicting the spatial patterns was found to be satisfactory. Both the NBN records and our model prediction suggest a wide distribution of ticks across mainland Scotland. Moreover, the model achieved ‘correct’ predictions (if the field NIP values fell within the range of simulated NIP values) in 22 forest sites (out of 24) sampled the field study [24], with two underestimations in Inverness (The Highlands). The lack of further field evidence on DIN or NIP, as well as the lack of sampling-level conditions of habitat and transmission host, prevented us from making further comparisons.

    3.1.3. Model sensitivity

    Sensitivity analysis of all model parameters was undertaken and compared with the results of previous model sensitivity investigations from the original model, which we extended in this study. These were found to be identical for all parameters [21,22]. In summary, the model had a relatively high sensitivity to tick mortality rate in the development phase from engorged larvae into questing nymphs, systemic transmission efficiencies, the mortality rate of questing nymphs, basal mortality rate of feeding larvae, transmission host feeding capacity for larvae, transmission host finding probabilities for larvae and nymphs, proportion of time-step spent in grassland for deer and the density of transmission hosts.

    3.2. Spatial pattern and seasonality of Lyme disease risk in Scotland

    The simulated spatial pattern of DIN in woodland is shown in figure 2a. The relative DIN, or the peak DIN (in week 37) as a percentage of the total DIN for all of Scotland, was used to analyse the capability of annual average temperature, elevation and deer density as simple indicators of high spatial disease risk. The relative DIN was predicted to decrease with increased elevation until it achieved a minimal value close to zero after 500 m (visual inspection of figure 2b; in accordance with empirical tick abundance data [51]). This was because, in the model, high altitude regions have (i) a colder environment which restricts tick development and questing activities (see [51]) and (ii) a large extent of non-woodland habitats with a low-transmission host density (e.g. in montane and inland rock), which was predicted to limit the chances of pathogen transmission (cf. adjacent habitat effects on pathogen transmission in Li et al. [21]). No consistent relationship was found between the predicted relative DIN and annual average temperature. However, where the annual average temperature was between 6°C and 10°C, the DIN was predicted to be lower in warmer places (visual inspection of figure 2c). Such places were predicted to be located in the central and southern parts of Scotland, and close to cities that had fewer deer, hence, the tick population was maintained at a low level. Otherwise (lower than 6°C or greater than 10°C), annual temperature did not predict the difference in DIN between cells, in which the effect of elevation was probably more dominant. Increasing deer density was predicted to be a good indicator of increasing disease risk in regions where deer density was between 15 and 22 heads km−2 (visual inspection of figure 2d). Alternatively, in places with deer density lower than 15 heads km−2 or greater than 22 heads km−2, the relative DIN was predicted to vary considerably. Some places had extremely high predicted DIN as they were suitable habitats for both deer and transmission hosts in low and medium altitude areas, while others had a predicted DIN of nearly zero as these cells were in high altitude areas, which are cooler and, therefore, have few ticks (after [51]) and fewer transmission hosts. The model predicted a similar seasonality of DIN for different Scottish regions, with some central and southern regions (Oban, Dumfries, Hamilton, Perth and Ayr) having a greater peak DIN value (figure 3a). A relatively higher risk was predicted in autumn (weeks 37–40) than in other seasons, lagging six weeks behind the peak weekly temperature (week 31–33, see the electronic supplementary material, figure S3.1f).

    Figure 2.

    Figure 2. Simulated spatial dynamics of Lyme disease risk in mainland Scotland. (a) Predicted spatial pattern of relative peak DIN (density of nymphs per square kilometre as per cent of total) and (bd) predicted influence of temperature increase on peak DIN on an increasing gradient of elevation (averaging interval of 1 m), annual average temperature (averaging interval of 0.1°C) and deer density (averaging interval of 0.1 heads km−2). (Online version in colour.)

    Figure 3.

    Figure 3. Predicted influence of temperature increase on the simulated seasonal variation of infected nymph abundance (summed DIN in woodland areas) in each AFRC (agriculture, food and rural communities) area office region. (Online version in colour.)

    3.3. Effects of climate warming on Lyme disease risk

    Regional climate warming led to a predicted expansion of the disease risk to higher altitudinal (figure 2b) and northern regions such as Thurso and Golspie (figure 3). With a temperature increase of up to 3°C, late winter and early spring were simulated as being no longer free of disease risk (figure 3). A warmer temperature was predicted to drive a greater proportion of ticks to resume development (from diapause) in spring and become more active over the winter. Thus, the duration of the interstadial development phase was predicted to be shortened and a greater tick population was predicted to survive before entering the questing phase. In addition, higher proportions of questing ticks were predicted to become active earlier in the spring and remain active later in the winter, resulting in a predicted prolonged duration of the tick questing season. Both could contribute to a greater frequency of tick–host contact and, hence, greater chance for pathogen transmission. By assuming that temperatures increased throughout the year, an upsurge in overall DIN was predicted for all regions (the electronic supplementary material, figure S1.1). Compared with the baseline pattern (no temperature increase), temperature increases of 1°C, 2°C and 3°C resulted in the predicted peak value of DIN increasing by 2, 7 and 11 times, respectively, and the extent of the predicted endemic area enlarging by 2.68%, 3.66% and 3.99%, respectively. The greatest increases in endemic area extent were predicted in Thainstone (3.96%, 6.24% and 7.80%), Hamilton (6.68%, 6.85% and 6.85%) and Ayr (5.12%, 5.56% and 5.60%). The overall effect of the temperature increases was predicted to be weaker as elevation increased (figure 2b). The predicted effect of climate warming scenarios was difficult to distinguish when plotted against annual average temperature (figure 2c) and deer density (figure 2d), because the two variables were both unable to explain the baseline DIN pattern (and the baseline DIN pattern was strongly correlated with the DIN increase in the climate warming scenarios; a significant Pearson's correlation coefficient of 0.67 was found between the predicted baseline DIN surface and the DIN surface under the temperature +1°C scenario).

    4. Discussion

    A mechanistic model was developed for mainland Scotland to predict the seasonal dynamics of the B. burgdorferi s.l. infected nymphal I. ricinus ticks (or DIN), a biophysical risk indicator of Lyme disease. The model integrates recent advances in simulating the temperature-driven temporal [1620] and spatial [21,22] dynamics of disease risk, with multi-disciplinary data based on the literature, empirical evidence, earth observation and model predictions. It was developed using an agent-based approach that consisted of three generalized grid layers (tick population, host population and landscape) and a range of transition rules describing their interactions under the influence of temperature. Considering the scale (the whole of mainland Scotland) and spatial resolution (1 km2) of this study, the model is more relevant for rodent-specialized Borrelia genospecies such as B. afzelii and B. burgdorferi s.s. than for the bird-specialized B. garinii or B. valaisiana because the densities of rodents are far greater than birds. However, recent empirical work in Scotland revealed that more that 40% of Borrelia is associated with birds, based on ticks collected across 25 woodland sites [3]. Further models specifically focusing on birds would be needed at a finer scale to gain a better picture of the Borrelia transmission network.

    It should be noted that the focus of this study was on long-term risk and the predicted DIN pattern, i.e. the results were based on 12-year (2000–2011) mean weekly temperature surfaces. Model predictions may have been different if different temperature data had been used. Furthermore, we focused on the relative values of DIN rather than on absolute values, which are difficult to measure in the field due to the low sampling efficiency of the flagging method which varies widely depending on the vegetation type and due to variation in questing tick behaviour related to weather conditions at the time of sampling. The predicted risk patterns suggest a spatial endemic foci in the Highlands (e.g. Oban and Inverness) and Tayside (e.g. Perth) and temporal foci in August and September, which are similar to the pattern of human incidence records [52,53]. In line with a recent field survey [51] in which ticks were found to be strongly and negatively associated with elevation on nine Scottish hills, our model simulations predict a general decrease in DIN with increased altitude in Scotland up to 500 m, above which the tick populations were close to zero. The ability of two other variables (annual mean temperature and deer density) to act as a proxy of disease risk was also analysed, but their explanatory capability may only be satisfactory within certain intervals: the DIN was predicted to be negatively associated with annual temperature between 6°C and 10°C and positively with deer density between 15 and 22 head km−2.

    Climate warming has been found to accelerate tick phenology [49] and contribute to their geographical expansion [54]. Our model predicted similar overall conclusions from applying scenarios of regional warming to the Lyme disease model. First, the model predicted an overall rise in DIN across the year with a peak in autumn. Climate warming was predicted to contribute to a greater frequency of tick–host contact and, therefore, greater chance for pathogen transmission. The consequent predicted DIN increases were greater in low altitude regions of Scotland where DIN had already been established at a relatively high level. Second, the model predicted that climate warming could increase the extent of tick infested areas, albeit marginally (less than 4% increase under all three climate warming scenarios) as the distribution of ticks in mainland Scotland is currently widespread [24]. Newly emerged DIN cells were predicted to follow an altitudinal direction, as warming turned cold high altitude places into areas suitable for maintaining an adequate population of questing ticks for a long enough period of the year. Most of these regions were inhabited by deer, which help to sustain tick populations.

    We have faced certain methodological and data challenges when applying such a complex modelling approach to a large area in the real world in line with those cited for complex socio-ecosystem modelling by Filatova et al. [55]. The most critical methodological challenge is how to represent the mechanisms of B. burgdorferi s.l. transmission as this involves multiple vector life stages, pathogen genospecies, host types and movement patterns. The potentially important mechanisms underpinning the interactions of system components are numerous and complex, but there is insufficient field evidence for their full quantification within our modelling framework. Hence, we focused on advancing the treatment of temperature effects within the model on tick population ecology and host population seasonality, while a number of simplifying assumptions were used for other parts of the system (e.g. generalizing pathogen species and host populations) to avoid introducing too many uncertainties and computation burdens. Data challenges were numerous, but, in general, data were available for model input and evaluation. We integrated data from diverse sources with different levels of detail, but it was difficult to assess any biases in the model inputs, particularly for transmission host distributions, which have been identified as an important disease risk component in previous modelling studies [3]. Detailed field data on the seasonality of both tick and pathogen populations were not available across a wide range of locations to undertake a thorough evaluation of model output. It is also unfortunate that the model outputs remain highly sensitive to the parameters for which the empirical data are the least reliable (e.g. systemic transmission efficiencies, transmission host finding probabilities). This data-poor condition has been noted as a common bottleneck in the modelling of vectors and diseases [56]. Therefore, the evaluation and parametrization experiments had to rely on the literature and discontinuous data from multiple sources. Further research is needed to fill these gaps either by improving data collection or by advancing empirical model predictions.

    This study focused on the biophysical component of the disease transmission system. Disease risk (density of infected nymphal ticks) was derived from those processes underlying vector survival and pathogen circulation in wildlife populations and, hence, could be considered as the ‘hazard’. However, a high ‘hazard’ does not always indicate a high human infection rate. Social factors shaping the pattern of land use for human outdoor activities, such as walking, forest ranging, hunting, hiking, scouting, orienteering and gardening [5759], are likely to highly influence the disease pattern [53,60]. Many of these outdoor activities are also likely to be influenced by climate change. Such activities expose humans to questing ticks and trigger pathogen spillover to the human population. To date, combined hazard and exposure risk models remain rare for vector-borne diseases, and, to our knowledge, have not been developed for tick-borne diseases. Future integration of exposure assessment may benefit from recent advances in land-use models, such as, agent-based socio-demographical models for estimating the spatial pattern of human trips to forests [61]. Application of such models to more plausible scenarios of future climate and socio-economic change would also constitute a notable extension of existing research. The simple temperature sensitivity analysis undertaken in this study demonstrated the importance of climate warming for disease risk, but more complex scenarios (such as those from regional climate models, where the temperature change varies both spatially and temporally, and those from land-use models, which include woodland expansion or contraction) are needed to better understand the pathways of disease risk and the uncertainties associated with these. Involving stakeholders in such future models and scenario development exercises (e.g. as undertaken in the CLIMSAVE project [62]) could significantly advance strategies for adaptive Lyme disease management under environmental change.

    Data accessibility

    The sources of datasets used in this modelling study have been presented in the electronic supplementary material (§S3).

    Authors' contributions

    S.L., P.A.H. and M.D.A.R. conceived and designed the study. L.G. contributed to the methods and data. S.L. performed the modelling experiments, interpreted the results and drafted the manuscript with contributions by all other authors. P.A.H., L.G. and M.D.A.R. edited the manuscript. All authors gave final approval for publication.

    Competing interests

    We declare we have no competing interests.

    Funding

    The research leading to these results has received funding from the European Community's Seventh Framework Programme (FP7/2007–2013) under grant agreement number 603416 (The IMPRESSIONS Project—Impacts and Risks from High-End Scenarios: Strategies for Innovative Solutions; www.impressions-project.eu).

    Acknowledgements

    The authors thank the two anonymous referees for their time and constructive comments on the manuscript.

    Footnotes

    Published by the Royal Society. All rights reserved.

    References