INTRODUCTION
Ticks are the most important vectors of zoonotic disease-causing pathogens in Europe, transmitting the tick-borne encephalitis (TBE) complex of viruses, Anaplasma phagocytophyllum, Babesia and Rickettsia species and Borrelia burgdorferi sensu lato, the complex of bacteria that cause Lyme borreliosis, amongst others. Ixodes ricinus L. ticks are particularly implicated in pathogen transmission because they are almost ubiquitous across Europe and are generalist feeders, which allows for pathogen transmission among different host species. Ixodes ricinus are increasing in number and range in many parts of Northern Europe [reviewed by Medlock et al. (Reference Medlock, Hansford, Bormane, Derdakova, Estrada-Peña, George, Golovljova, Jaenson, Jensen, Jensen, Kazimirova, Oteo, Papa, Pfister, Plantard, Randolph, Rizzoli, Santos-Silva, Sprong, Vial, Hendrickx, Zeller and Van Bortel2013)].
In any given geographical region tick population dynamics are dependent on a number of biotic and abiotic factors including the density of different host species, and other factors that influence survival and activity such as temperature and humidity and vegetation types, the latter of which provide habitats for different hosts and create different microclimates.
Mathematical models have been used extensively to predict the dynamics of tick populations under different conditions including climate change. However, high tick densities do not necessarily mean high prevalence or risk of tick-borne pathogens, since this is dependent not only ticks but also competent transmission hosts. Therefore, models have also been used to predict the tick-borne pathogen dynamics and the theoretical effectiveness of different tick-borne pathogen control methods under different environmental or management scenarios. In this paper, we will review the use of those models for different systems, summarize the key results in different contexts and discuss possible future directions of mathematical modelling of tick-borne pathogens.
MATHEMATICAL MODELS OF TICK POPULATION DYNAMICS
Although there are a number of different tick species globally this review will focus on I. ricinus and we will specify when we cite any papers which refer to other species.
The I. ricinus life cycle develops from the egg, through two immature stages (larvae and nymph) to the adult stage. Each immature stage requires a blood meal from a suitable vertebrate host before developing to the next stage and the adult female requires a blood meal before producing eggs. Adult females feed primarily on large mammals such as deer, sheep or hares whilst the immature stages can also feed on smaller vertebrates such as mice, voles and birds (e.g. Gray, Reference Gray1998). The I. ricinus life cycle usually takes 3–4 years to complete (Fig. 1).
Fig. 1. Schematic diagram of the I. ricinus life cycle with the type of host that they are able to feed on at each stage.
In winter ticks often enter behavioural diapause induced by cold and/or short day length (Randolph et al. Reference Randolph, Green, Hoodless and Peacey2002; but see Gray, Reference Gray1987). Therefore, tick activity is highly seasonal with ticks in Northern Europe being active mainly between spring and autumn when temperatures are warm enough. Activity is inhibited by cold temperatures but increases with temperature up to a limit [12–20 °C depending on population e.g. Gilbert et al. (Reference Gilbert, Aungier and Tomkins2014); Tomkins et al. (Reference Tomkins, Aungier, Hazel and Gilbert2014)]. Tick host-seeking (questing) activity can also be inhibited by low relative humidity or high saturation deficits [this is a function of relative humidity and temperature and gives an estimate of the drying power of the air; Perret et al. (Reference Perret, Guigoz, Rais and Gern2000)]. After feeding, ticks also become inactive due to physiological diapause while they develop into the next stage (Randolph et al. Reference Randolph, Green, Hoodless and Peacey2002).
One of the first mathematical models developed to describe tick population dynamics was published in 1981 (Gardiner et al. Reference Gardiner, Gettinby and Gray1981). This study used empirical data from experiments to predict how tick development times depend on temperature. They did not put this into a formal predictive modelling framework but they did try to determine functional relationships between development time and different measures of temperature (i.e. air and soil temperature). In particular, they looked at how experimentally predicted development times estimated in the laboratory translated to the field where temperature fluctuations are much less predictable. They found that soil temperatures recorded at a depth of 50 mm are useful predictors for larval and nymphal development phases. In terms of egg development time they found that air temperatures are useful for predicting the development time of eggs laid in the spring but soil temperature is a better predictor for those laid in autumn. They suggested that this might be because during diapause eggs may be conditioned to develop according to the temperature of their environment rather than air temperature.
Mount and Haile (Reference Mount and Haile1989) developed a computer simulation model of the American dog tick Dermacentor variabilis (Say). This model simulated the effects of environmental variables such as ambient temperature, habitat and host density on American dog tick population dynamics. They validated the model by comparing its predictions with empirical data from Virginia, Maryland and Massachusetts. The authors concluded that the model produced acceptable values for equilibrium population densities and seasonal activity patterns and went on to extend this model to include Rocky Mountain spotted fever dynamics (Cooksey et al. Reference Cooksey, Haile and Mount1990).
Over the last 40 years Sarah Randolph and collaborators have written a large number of papers on tick biology and population dynamics. These are largely empirical; however, there are also some which model tick population dynamics. The first of these came in 1997 (Randolph and Rogers, Reference Randolph and Rogers1997) where they presented a simulation model of the African tick Rhipicephalus appendiculatus. This simulation model incorporated temperature-dependent rates of egg production and development, climate-driven density-independent mortality rates and density-dependent regulation of both nymphs and adults. The model successfully described both the seasonality and annual range of variation in numbers of each tick stage observed at each of four test sites in Uganda, Burundi and South Africa.
In 2002, Randolph et al. (Reference Randolph, Green, Hoodless and Peacey2002) used empirical data on tick counts, various microclimatic factors and fat contents of ticks to create a population model explaining seasonality of I. ricinus in the UK. This study showed large variation in questing activity between years, but the date of questing (i.e. host-seeking activity) in 1 year was used to predict the start of questing for the next stage the following year, with reasonable accuracy. This was an important paper that also found evidence of two cohorts of ticks within a life stage within a season. Those nymphs with higher relative fat contents had emerged and become active more recently than those with lower fat contents. The suggestion was that spring-questing nymphs had overwintered, having fed as larvae the previous late summer or autumn; meanwhile autumn-questing nymphs had fed as larvae in the spring of the same calendar year.
More recently, Dobson et al. (Reference Dobson, Finnie and Randolph2011b ) used a stage-classified Leslie matrix model to break the tick life cycle into the key parts, with a particular focus on two types of diapause: developmental and behavioural, with the latter being important in determining how many times a year an individual tick might feed. This model was then used by Dobson and Randolph (Reference Dobson and Randolph2011a ) to make long-term predictions of the effects of host densities, climate and acaricide treatment of hosts on tick populations.
In 2005, Ogden et al. (Reference Ogden, Bigras-Poulin, O'Callaghan, Barker, Lindsay, Maarouf, Smoyer-omic, Waltner-Toews and Charron2005) developed a model of Ixodes scapularis Say (1821) in which tick development rates were modelled as temperature-dependent time delays. Time spent in egg and engorged tick states and questing activities were all temperature dependent. The parameters were estimated using data taken from Ogden et al. (Reference Ogden, Lindsay, Charron, Beauchamp, Maarouf, O'Callaghan, Waltner-Tiews and Barker2004). The model was validated using data from Ontario and Maryland and in both cases the observed seasonal activity patterns were predicted by the model. The models were then used to predict theoretical geographical limits for the establishment of I. scapularis in Canada. The model predicted that the temperature conditions which are suitable for the tick are wider than the existing distribution, implying that there is potential for spread.
An age-structured stochastic model was used to describe the dynamics of tick populations by Hancock et al. (Reference Hancock, Brackley and Palmer2011). They focused on the effect of temperature on the development between each stage of the tick life cycle, i.e. from egg to larva, larva to nymph, nymph to adult and adult laying eggs. This model also introduced pathogen dynamics into the model. This allowed the model to predict that, if a pathogen is introduced into the system, it is most likely to persist if it is introduced at a time of year of peak tick questing.
A completely different approach was adopted by Schwarz et al. (Reference Schwarz, Maier, Kistemann and Kampen2009) who used statistical methods to identify the relationship between vegetation and tick distribution. Ixodes ricinus tick count data were correlated with plant communities, and the resulting relationship was used to predict I. ricinus distribution across the German nature reserve Siebengebirge, using Geographic Information Systems (GIS). A similar process was undertaken by Braga (Reference Braga2012) to identify the associations between habitat, host densities, temperature and other climatic factors on observed tick abundance at sites across Scotland. The resulting output was used to predict tick abundance over all of Scotland according to GIS-based environmental data, and visualized as a series of raster maps showing predicted tick abundance. The key parameters in this basic algorithm were then altered in accordance with environmental change projections (climate change and woodland expansion), to produce predictions of future tick abundance over Scotland due to environmental change scenarios.
Jore et al. (Reference Jore, Viljugrein, Hofshagen, Brun-Hansen, Kristoffersen, Nygård, Brun, Ottesen, Sævik and Ytrehus2011) also used a statistical method to investigate I. ricinus tick dynamics. A principle component analysis provided a model which explained 67% of the variation in past I. ricinus densities in Norway. The study suggests that I. ricinus have expanded northwards since 1983.
Summary
For almost 35 years mathematical models of tick dynamics have been developed. The models have largely focused on the impact of environmental factors on these dynamics. Field observations show that tick life stages emerge at different points in the season and peak at different times in different geographical regions. In some areas, we can have bimodal tick dynamics within a year (e.g. Tagliapietra et al. Reference Tagliapietra, Rosa, Arnoldi, Cagnacci, Capelli, Montarsi, Hauffe and Rizzoli2011) and in other areas there is only one peak. The models described above have been able to replicate the observed tick dynamics for particular geographical areas, tick species and environmental conditions. However, it is clear that in order to be able to predict tick dynamics we would need to have key pieces of information about the environment (and particularly the temperature) in which they live.
Lorenz et al. (Reference Lorenz, Dhingra, Chang, Bisanzio, Liu and Remais2014) explicitly looked at the extrapolation of landscape model results based on one area or time to other spatial or temporal systems for Lyme disease and I. scapularis and concluded that models based on measures of vegetation, habitat patch characteristics and herbaceous landcover emerged as effective predictors of observed disease and vector distribution. These would therefore be important characteristics of an area to measure in order to predict these distributions.
MATHEMATICAL MODELS OF TICK-BORNE PATHOGEN DYNAMICS
Modelling of tick-borne pathogens has focussed on a small number of pathogens which are important for human or animal health and welfare. The three main systems which have been modelled extensively are Louping ill virus (LIV), western tick-borne encephalitis virus (TBEV) and B. burgdorferi sensu lato, the causative agent of Lyme disease. This section will focus largely on LIV since this pathogen has the largest body of modelling work and it is the area of expertise of the authors. It also illustrates many of the biological features which need to be incorporated into models and so is a good case study for models of other system.
In general, transmission of these pathogens can occur in three ways [although also see Park et al. (Reference Park, Robertson, Campbell, Foster, Russell, Newborn and Hudson2001) discussed below for LIV]. The most common form of transmission occurs when susceptible ticks feed on infected hosts who have virus in their bloodstream (viraemic hosts) and pick up the virus. These ticks then moult into their next developmental stage and when they take their next blood meal then they can pass the pathogen onto a susceptible host, this will be a different individual and can also be a different host species (Labuda and Nuttall, Reference Labuda and Nuttall2004). The second method is vertical transmission; for some pathogens infection is passed from adult ticks to eggs and onto larvae (Labuda and Nuttall, Reference Labuda and Nuttall2004). Finally, for some hosts and some pathogens there can be non-viraemic or co-feeding transmission in which susceptible ticks feeding near to infectious ticks can pick up infection without the host having a viraemic response (Jones et al. Reference Jones, Davies, Steele and Nuttall1987).
Louping ill virus
A large number of increasingly complex models have been used to help us understand LIV, which is the western-most variant of Western TBEV. LIV is transmitted by I. ricinus and causes disease in livestock, especially sheep Ovus aries, as well as red grouse, Lagopus lagopus scoticus, a valuable game bird. A vaccine has been developed for livestock but not for red grouse that are highly susceptible to the disease, with 78% mortality rates in experimentally infected birds in the laboratory (Reid, Reference Reid and Wilde1976). The hosts and transmission cycle of this complex virus system has been recently reviewed (Gilbert, Reference Gilbert2015), but mathematical models can be extremely useful in helping to identify gaps in our biological knowledge of the system, identifying the relative importance of different host species hosts, and predicting the effectiveness of potential control strategies.
The first mathematical model of LIV was presented by Hudson et al. (Reference Hudson, Norman, Laurenson, Newborn, Gaunt, Gould, Reid, Bowers and Dobson1995), where a series of coupled ordinary differential equations describing LIV on red grouse moorland was presented. This model explored the interactions between ticks and red grouse and their role in the dynamics of LIV. The model predicted that grouse alone cannot support a tick population since very few adult ticks feed on grouse, therefore other hosts are required to complete the tick life cycle. Within this model the alternative hosts were mountain hares Lepus timidus, although similar later studies examined the role of red deer Cervus elaphus (Gilbert et al. Reference Gilbert, Norman, Laurenson, Reid and Hudson2001; Norman et al. Reference Norman, Ross, Laurenson and Hudson2004) and sheep (Porter et al. Reference Porter, Norman and Gilbert2011). Hudson et al. (Reference Hudson, Norman, Laurenson, Newborn, Gaunt, Gould, Reid, Bowers and Dobson1995) also calculated a formula for the conditions for persistence of both ticks and LIV. For tick persistence a sufficient number of hosts (or combination of host types) which can feed all stages of ticks are required, while LIV persistence also requires a competent LIV transmission host (red grouse in this model) to make up a sufficient proportion of the total tick hosts. This means that, in order for the pathogen to persist one needs enough tick hosts to maintain the tick population, with a sufficient number of these being pathogen-transmitting hosts. This threshold formula comes from the basic reproductive rate or number, R0, which is a concept used widely in epidemic models. When R0 > 1 then the pathogen persists and when R0 < 1 the pathogen dies out.
Some more complex later LIV models have also predicted an eventual ‘dilution effect’ where pathogen prevalence declines if there are too many non-pathogen-transmitting tick hosts (hosts which do not transmit the pathogen such as deer) compared with competent transmission hosts which causes potential pathogen-transmitting bites to be ‘wasted’ and the effect of the pathogen to be diluted (Norman et al. Reference Norman, Bowers, Begon and Hudson1999; Gilbert et al. Reference Gilbert, Norman, Laurenson, Reid and Hudson2001).
Sheep are known to produce a LIV viraemia after infection, and are known to be competent transmission hosts. However, the role of lambs is less well understood; if ewes have been bitten by infected ticks, their young lambs acquire immunity from the virus from drinking the colostrum from their mothers in the first few days or weeks of life. However, as the lambs age this immunity wanes, leaving them at risk of contracting LIV. Thus, lambs could potentially have a role as a reservoir host. Therefore, another differential equation model was created to understand the role that lambs may play as a reservoir of LIV. The model predicted that, whilst in theory large numbers of lambs could act as a reservoir for the virus, it is more likely that, in most situations, these numbers are probably small (Laurenson et al. Reference Laurenson, Norman, Reid, Pow, Newborn and Hudson2000).
Laurenson et al. (Reference Laurenson, Norman, Gilbert, Reid and Hudson2003) examined the impact of near-eradication of mountain hares on tick burdens and LIV seroprevalence in red grouse, using both empirical data and differential equation models. The models compared the scenario where mountain hares simply act as tick amplifying hosts to a scenario where hares were both tick hosts and non-viraemic transmission hosts. It was found that the model which included non-viraemic transmission produced predictions that fitted the data better than the simpler model did. Laboratory experiments had already identified mountain hares as competent transmission hosts (through supporting non-viraemic transmission between co-feeding ticks) in the laboratory (Nuttall and Jones, Reference Nuttall, Jones, Dusbabek and Buvka1991; Jones et al. Reference Jones, Gaunt, Hails, Laurenson, Hudson, Reid, Henbest and Gould1997). In addition, models have shown that non-viraemic transmission via co-feeding may allow the virus to persist more readily than it would otherwise have done, and allow the virus to persist even in the absence of viraemic hosts if the level of non-viraemic transmission is high enough (Norman et al. Reference Norman, Ross, Laurenson and Hudson2004). However, the Laurenson et al. (Reference Laurenson, Norman, Gilbert, Reid and Hudson2003) study was important in demonstrating that mountain hares can be LIV reservoir hosts in the field. There were large management repercussions to this research, as many grouse moor managers over Scotland began large-scale culls of mountain hares, leading to political issues (reviewed by Harrison et al. Reference Harrison, Newey, Gilbert, Haydon and Thirgood2010; Gilbert Reference Gilbert2015). Models again had political impact by providing evidence against culling mountain hares: while the Laurenson et al. (Reference Laurenson, Norman, Gilbert, Reid and Hudson2003) system included only red grouse and mountain hares, most areas in Scotland managed for grouse hunting also have deer. Therefore, Gilbert et al. (Reference Gilbert, Norman, Laurenson, Reid and Hudson2001) modelled a three-host system, including deer as well as red grouse and mountain hares. Importantly, this three-host model predicted that LIV would always persist in the presence of even low densities of deer, even if all mountain hares were culled. This was because red grouse are transmission hosts for the virus while deer, although not competent transmission hosts, are important hosts for all stages of tick, so together both virus and tick life cycles can be maintained. This Gilbert et al. (Reference Gilbert, Norman, Laurenson, Reid and Hudson2001) model has been crucial in the arguments against large-scale mountain hare culls (Harrison et al. Reference Harrison, Newey, Gilbert, Haydon and Thirgood2010; Gilbert, Reference Gilbert2015).
Mathematical models have also been used in helping identify which pathogen control methods could be theoretically most effective in LIV control. Porter et al. (Reference Porter, Norman and Gilbert2011) developed models to predict the effectiveness of using acaricide-treated sheep as a tool to control ticks and LIV in red grouse. The model predicted that the presence of deer limits the effectiveness of such a strategy, but for certain conditions the use of acaricide on sheep could theoretically be a viable method for controlling ticks and LIV providing that high numbers of sheep are treated and acaricide efficacy remains high, while deer densities must be very low (Porter et al. Reference Porter, Norman and Gilbert2011). Due to this predicted adverse impact of deer on the success of treating sheep to control ticks and LIV, and because deer are known to maintain high tick population densities in Scotland and move ticks between habitats (Ruiz-Fons and Gilbert, Reference Ruiz-Fons and Gilbert2010; Jones et al. Reference Jones, Webb, Ruiz-Fons, Albon and Gilbert2011; Gilbert et al. Reference Gilbert, Maffey, Ramsay and Hester2012), models were then developed to test the theoretical effectiveness of acaricide-treated deer on controlling ticks and LIV (Porter et al. Reference Porter, Norman and Gilbert2013a ). The model predicted that treating deer could control ticks and LIV if high acaricide efficacies were maintained and if a large proportion of the deer population was treated. Furthermore, effectiveness was improved if there were only low densities of deer. However, although the model predicted that this control method is theoretically plausible, it is unlikely that the conditions could be met in practical terms, in wild deer. Therefore, using an age-structured differential equation model, including splitting the grouse life cycle to represent the different behaviour between chicks and adults, Porter et al. (Reference Porter, Norman and Gilbert2013b ) investigated whether acaricide treatment of the grouse themselves could help reduce ticks in the environment and LIV in the grouse population. Again, this was theoretically possible, but in the presence of deer, high acaricide efficacies were required and high proportions of the grouse population treated, were needed for successful control. This is due to the deer amplifying the tick population. These types of models can therefore be of use in decision-making by land managers for choosing disease control options, such as whether to try a certain control method or not depending on the situation in a specific area, taking into account any practical difficulties.
It is generally assumed that LIV is transmitted through ticks biting their hosts, and model parameterization generally reflects this assumption. However, red grouse chicks frequently eat invertebrates, including ticks (Park et al. Reference Park, Robertson, Campbell, Foster, Russell, Newborn and Hudson2001). This is a potentially important route of transmission: it has been suggested that 73–98% of LIV infection in red grouse in their 5 year could stem from ingestion (Gilbert et al. Reference Gilbert, Jones, Laurenson, Gould, Reid and Hudson2004). Introducing this infection route to LIV modelling has an interesting effect: when using the standard method for calculating the basic reproduction number for the persistence of LIV, then the algebraic results and numerical simulations do not match. The standard method of analysis causes virus persistence to be underestimated, as the ingestion of infected ticks causes a feedback loop where the virus can persist with seemingly insufficient hosts (Porter et al. Reference Porter, Norman and Gilbert2011). This phenomenon requires further investigation, as it may indicate interesting gaps in our knowledge of the biology of the LIV system as well as an anomaly in the current modelling approach.
In the LIV models described above, there has been no explicit spatial component to the models. However, Watts et al. (Reference Watts, Palmer, Bowman, Irvine, Smith and Travis2009) investigated the interaction between neighbouring areas by expanding the previously existing LIV models into a two-patch system with host movement between patches. Comparison with empirical data showed that whilst the one-patch model was a reasonable indicator for tick numbers, it tended to underestimate the prevalence of the LIV. When considering the two-patch model, the results depended largely on finding the appropriate balance of deer movement between the two sites (Watts et al. Reference Watts, Palmer, Bowman, Irvine, Smith and Travis2009). Jones et al. (Reference Jones, Webb, Ruiz-Fons, Albon and Gilbert2011) developed a different type of differential equation model, which explicitly tracked the number of ticks on each host, to predict how deer moving ticks from forest onto moorland might affect ticks and LIV in red grouse on the moorland. The assumption was that ticks are more abundant in forest than on moorland, which is supported by empirical data (Ruiz-Fons and Gilbert, Reference Ruiz-Fons and Gilbert2010). This model predicted the highest levels of LIV in moorland to occur where it is bordering forest regions, due to higher tick numbers there. Furthermore, this model was important in examining for the first time the impact of landscape heterogeneity on predicted pathogen levels: virus prevalence was predicted to be higher in landscapes that have larger forest patches, and higher landscape fragmentation, which increases the number of borders between the two habitats (Jones et al. Reference Jones, Webb, Ruiz-Fons, Albon and Gilbert2011).
Summary
The transmission, persistence and dynamics of LIV are complex with many interacting factors to take into account. The focus of the modelling work described above has been on trying to understand the roles that different hosts play in maintaining these dynamics. Hosts can play three possible roles, they can either simply act as tick amplifiers (e.g. deer) or they can both amplify ticks and transmit virus (e.g. sheep for viraemic transmission or hares for non-viraemic transmission) or finally they can transmit the disease but not support the ticks (e.g. grouse). The ability to control the virus in any particular geographical area is highly dependent on the host community structure. In addition there are practical issues involved in trying to control the virus in this system which is made up of mostly wild hosts. There are both physical difficulties in delivering treatment and legislative difficulties in which treatments are permitted.
Other tick-borne pathogens
Tick-borne encephalitis
TBE is a neurological disease which is of significant public health interest across mainland Europe. It is caused by the TBEV, which is primarily transmitted by I. ricinus ticks, where rodents act as the competent host for the virus.
There are two significant ways in which deer can influence TBEV dynamics. Firstly, as deer are the main host which I. ricinus adults feed on, their presence, as with LIV, has an amplification effect on tick abundance. Secondly, as deer do not support TBEV transmission, very high deer densities can eventually lead to the dilution effect lowering TBEV levels (again similar to model predictions of LIV).
In both 2003 and 2007, Rosa and co-authors extended the models of Norman et al. (Reference Norman, Bowers, Begon and Hudson1999) to explicitly include the questing and feeding tick stages and the aggregation of ticks on the hosts. They investigated changes in host densities and different infection pathways to determine when the dilution effect might occur. They found the new result that the dilution effect might occur at high densities of disease competent hosts. The authors state that better information on tick demography would be needed before it would be possible to predict whether this effect would happen in the field. However, there is some evidence that this is the case in the TBE system (Perkins, Reference Perkins2003).
In 2012, the same Italian group published a pair of papers taking both an empirical and theoretical approach to understanding the effect of deer density on tick distributions on rodents and therefore on the risk of TBE in a region. Cagnacci et al. (Reference Cagnacci, Bolzoni, Rosa, Carpi, Hauffe, Valent, Tagliapietra, Kazimirova, Koci, Stanko, Lukan, Henttonen and Rizzoli2012) empirically found a hump-shaped relationship between deer density and ticks feeding on rodents, and a negative relationship between deer density and TBE occurrence. Twinned with this, a model was developed by Bolzoni et al. (Reference Bolzoni, Rosa, Cagnacci and Rizzoli2012) to explain these findings. They found hump-shaped relationships between deer density and both the number of ticks feeding on rodents and TBEV prevalence in ticks. For low deer densities this can be explained by the tick amplification effect, for high deer densities the virus dilution mechanism dominates the dynamics.
The role of climate change on tick-borne pathogen prevalence was scrutinized by Randolph (Reference Randolph2008). In this study, TBEV was used as a case example. A statistical model was used to show that climate change is not enough to explain historical changes in TBE incidence within Europe. An alternative model was presented, showing how the introduction of further factors allowed for a better model fit of the data. Crucially, such a model included socioeconomic factors such as unemployment, agricultural practices and income. Zeman et al. (Reference Zeman, Pazdiora and Benes2010) used GIS analysis to similarly find that heterogeneity in TBE trends cannot be fully explained by geographic and climatic factors. However, they also found that the inclusion of socioeconomic conditions could not satisfactorily explain the anomalies.
Summary
As with Louping ill the persistence and dynamics of TBE are dependent on host densities and deer play a crucial role in this. Some of the papers described above, particularly in the 2003 and 2007 papers, Rosa et al. (Reference Rosa, Pugliese, Norman and Hudson2003, Reference Rosa and Pugliese2007) present general results which could apply to a number of different tick-borne pathogens and, in particular the results that dilution effects are very dependent on tick demography and density-dependent constraints are true more generally than just for TBE. In most of the models presented here TBE has been a case study of a model which addresses more general questions.
Lyme disease
Borrelia burgdorferi s.l. is the suite of spirochete bacteria which causes Lyme disease. This is a pathogen which has a wildlife reservoir but infects humans in the northern hemisphere.
Porco (Reference Porco1999) used a time-independent differential equation model to investigate how the prevalence of B. burgdorferi s.l. in I. scapularis (Say) nymphs is affected by various model parameters. The infectivity of white-footed mice Peromyscus leucopus (a key transmission host in the Eastern USA) was predicted to be the parameter which increased B. burgdorferi s.l. prevalence the most, whilst a tenfold increase in the density of deer (which do not transmit the pathogen) significantly reduced B. burgdorferi s.l. prevalence, suggesting that this is another system where the dilution effect can occur.
Zhang and Zhao (Reference Zhang and Zhao2013) presented a seasonal reaction–diffusion model of Lyme disease, utilizing it to study the dynamics of the system in bounded and unbounded spaces. For bounded habitats a threshold for pathogen persistence was predicted, whilst for unbounded habitats they were able to predict the speed of pathogen spread.
In their 2007 paper, Ogden et al. (Reference Ogden, Bigras-Poulin, O'Callaghan, Barker, Kurtenbach, Lindsay and Charron2007) considered the work of Wilson and Spielman (Reference Wilson and Spielman1985) and hypothesized that the transmission cycles of B. burgdorferi are very efficient in north-eastern North America because the seasonal activity of nymphal and larval I. scapularis is asynchronous. They then developed a simulation model which integrated transmission patterns imposed by seasonal asynchronous nymph and larvae with a model of infection in white footed mice. They parameterized the model for B. burgdorferi and Anaplasma phagocytophilum as examples. They found that duration of host infectivity, transmission efficiency to ticks and co-feeding transmission are the major factors determining fitness of pathogens in I. scapularis in North America.
The same group then wrote a series of papers looking I. scapularis in Canada where is established in some places and emerging in others. In Wu et al. (Reference Wu, Duvvuri, Lou, Ogden, Pelcat and Wu2013) they developed a temperature driven map of the basic reproductive number for the ticks and found that while the geographical extent of suitable tick habitat is expected to increase with climate warming the rate of invasion will also increase. In a subsequent paper, Ogden et al. (Reference Ogden, Lindsay and Leighton2013) investigated the speed of B. burgdorferi invasion after establishment of ticks. The model showed that the number of immigrating ticks was a key determinant of pathogen invasion and so the authors hypothesized that a 5-year gap would occur between tick and B. burgdorferi invasion in Eastern Canada but a much shorter gap in Central Canada. This was consistent with empirical evidence.
Summary
Borrelia burgdorferi is another pathogen for which the dilution effect appears to occur. In this case rodents are the main reservoir host and B. burgdorferi is emerging in a number of different areas as the tick hosts expand their range in response to climate change or socioeconomic factors.
More general models of tick-borne pathogen
More generally, Hartemink et al. (Reference Hartemink, Randolph, Davis and Heesterbeek2008) determined ways of characterizing the basic reproductive number in a tick-borne pathogen system which has multiple transmission routes using the next generation matrix (e.g. Diekmann et al. Reference Diekmann, Heesterbeek and Roberts2010). They showed that the complexities of the tick transmission cycle can be overcome by separating the host population into epidemiologically different types of individuals and constructing a matrix of reproduction numbers. They then used field and experimental data to parameterize this next-generation matrix for B. burgdorferi s.l. and TBEV.
Dunn et al. (Reference Dunn, Davis, Staecy and Diuk-Wasser2013) used a general model of tick-borne pathogens to study the basic reproductive number and found that the transmission efficiency to the ticks, the survival rate from feeding larvae to feeding nymphs and the fraction of nypmhs to find a competent host are the most important factors in determining R0.
Another general tick-borne pathogen model was created by Zeman (Reference Zeman1997), where reported cases of disease were smoothed over to create risk maps for Lyme disease and TBE in Central Bavaria. This study indicated that B. burgdorferi s.l. is wider spread than TBEV, but that both pathogens share the same main foci. Similarly, Hönig et al. (Reference Hönig, Švec, Masař and Grubhoffer2011) assessed the suitability of various habitats for supporting I. ricinus ticks, creating a model with which they were able to create a tick-borne pathogen risk map for South Bohemia, which was compared to clinical cases of TBE for validation. The model suggested that the areas mostly suitable for tick-borne pathogens were along the river valleys. However, when human activity is taken into account, the surroundings of large settlements are equally likely to provide tick-borne pathogen cases.
Another aspect of transmission which is considerably less well understood is the pattern of aggregation of ticks on hosts. Ferreri et al. (Reference Ferreri, Giacobini, bajardi, Bertolotti, Bolzoni, Tagliapietre, Rizzoli and Rosa2014) analysed a 9-year time series of I. ricinus feeding on Apodemus flavicollis mice, the reservoir host for TBE in Trentino, Northern Italy. The tail of the distribution of the number of ticks per host was fitted to three theoretical distributions. The impact of these distributions on pathogen transmission was investigated using a stochastic model. Model simulations showed that there were different outcomes of disease spread with different distribution laws amongst ticks and so it is important to have data on these distributions which can vary seasonally.
The models discussed above are not an exhaustive list, but do describe models which help us to understand many of the different complexities of tick-borne pathogen systems, and showcase the diversity of models now being developed for a wide range of end uses.
KNOWLEDGE GAPS AND FUTURE DIRECTIONS
As we have seen mathematical models have been used for more than 30 years to help to predict tick dynamics and subsequently pathogen dynamics. The models presented here have been used in two ways, firstly to predict when tick densities are at their peak within a year and how that peak varies with environmental factors. Secondly, they have been used to predict pathogen persistence for different combinations of available host species with different transmission competencies. In particular, they have looked at the interaction between tick amplifying hosts and disease-transmitting hosts and how densities of these hosts could be manipulated to control the disease.
One of the problems of these modelling studies is the difficulty in gathering empirical data to validate the model results. This is largely because there is a great deal of variability between sites in terms of habitat cover, microclimate and host densities. This is not unique to the tick system; it is difficult for a number of reasons to carry out experiments in natural systems. It is also difficult to measure realistic tick densities (e.g. Dobson, Reference Dobson2014).
However, most of the models described here have succeeded in doing some type of validation and they provide useful qualitative results.
Future modelling approaches are likely to be focussed in two main areas. One is to look at spatial patterns of tick and disease risk, and in particular to link environmental information in GIS systems to models of tick and pathogen dynamics in a mechanistic way. These models can then be used to predict the impact of climate change on tick and disease risk across a given geographical region. This type of modelling is currently being carried out by the authors for both LIV and Lyme disease in Scotland. The advantage of this type of modelling is that it is generalizable and could be applied to any country with the right type of environmental data available in GIS form. It can also predict risks are going to change over time rather than only looking at the end points as has been done before (e.g. Braga, Reference Braga2012).
If we can identify which areas are going to have significant increases in disease risk, then we can inform policy makers and target control efforts. For example, if we could identify which areas are going to have higher and lower Lyme disease risk then we could target efforts to educate the public on how to avoid being bitten in those high risk areas.
The second direction in which we predict tick modelling will move is to further develop a new modelling technique which was introduced in Jones et al. (Reference Jones, Webb, Ruiz-Fons, Albon and Gilbert2011). In that paper, the authors developed a model which keeps track of the number of hosts with a particular number of ticks on it. This allows us to consider the distribution of ticks amongst the hosts and the spatial element of the model allowed hosts to be infected in one place and for ticks to drop off them in another area. This allows us to model the movement of ticks which is not possible in the more traditional models. This ability adds to the level of biological realism and allows us to answer questions which it is not normally possible to answer. For example, we would be able to model control more explicitly in terms of how many ticks are lost from an individual host. This approach allows us to relax the intrinsic assumption of most models that all host individuals, and indeed, all tick individual behave in the same average way and to look at distributions of behaviour which should bring further insights into these important systems.
Whilst it is interesting to speculate about future directions in the modelling of tick-borne diseases, it is also likely that future modelling will be developed in response to new insights, either in the current systems, or in emerging diseases which become a problem in the light of climate change, for example. New models or new methods of analysis will have to be developed in order to answer specific questions for particular biological systems.
ACKNOWLEDGEMENT
AJW is supported by a Partnership Impact Ph.D. Studentship funded by both the University of Stirling and the James Hutton Institute.
