Pattern of tick aggregation on mice: larger than expected distribution tail enhances the spread of tick-borne pathogens

P attern of tic k aggregation on mice: larger than exp ected distribution tail enhances the spread of tic k-b orne pathogens Luca F erreri 1 , 2 , ∗ , Mario Giacobini 1 , 2 , 3 , P aolo Ba jardi 1 , 2 , Luigi Bertolotti 1 , Luca Bolzoni 4 , 5 , V alen tina T agliapietra 5 , Annapaola Rizzoli 5 , Rob erto Ros` a 5 1 Computational Epidemiology Group, Departmen t of V eterinary Sciences, Univ ersit y of T orino, Italy 2 Applied Researc h on Computational Complex Systems Group, Departmen t of Computer Science, Universit y of T orino, Italy 3 Complex Systems Unit, Molecular Biotec hnology Centre, Univ ersit y of T orino, Italy 4 Istituto Zo oprofilattico Sp erimen tale della Lom bardia e dell’Emilia Romagna, P arma, Italy 5 Dipartimen to Bio div ersit` a ed Ecologia Molecolare, Cen tro Ricerca e Inno v azione, F ondazione Edm und Mac h, San Michele all’Adige (TN), Italy ∗ E-mail: luca.ferreri@unito.it Abstract The spread of tic k-borne pathogens represen ts an imp ortan t threat to human and animal health in many parts of Eurasia. Here, w e analysed a 9-y ear time series of Ixo des ricinus ticks feeding on Ap o demus flavic ol lis mice (main reserv oir-comp eten t host for tic k-b orne encephalitis, TBE) sampled in T rentino (Northern Italy). The tail of the distribution of the num b er of tic ks p er host was fitted b y three theoretical distributions: Negativ e Binomial (NB), Poisson-LogNormal (PoiLN), and P ow er-La w (PL). The fit with theoretical distributions indicated that the tail of the tick infestation pattern on mice is b etter describ ed b y the PL distribution. Moreo ver, we found that the tail of the distribution significantly c hanges with seasonal v ariations in host abundance. In order to inv estigate the effect of differen t tails of tic k distribution on the in v asion of a non-systemically transmitted pathogen, we sim ulated the transmission of a TBE-like virus b et ween susceptible and infective tic ks using a stochastic mo del. Mo del simulations indicated different outcomes of disease spreading when considering differen t distribution la ws of tic ks among hosts. Sp ecifically , we found that the epidemic threshold and the prev alence equilibria obtained in epidemiological simulations with PL distribution are a go od appro ximation of those observ ed in simulations feed by the empirical distribution. Moreo ver, w e also found that the epidemic threshold for disease in v asion w as lo wer when considering the seasonal v ariation of tick aggregation. Autho r Summary Our w ork analyses a 9-year time series of tick co-feeding patterns on Y ellow-nec k ed mice. Our data shows a strong heterogeneit y , where most mice are parasitised by a small n um- b er of ticks while few host a m uch larger num b er. W e describ e the num b er of ticks 1 p er host by the commonly used Negative Binomial mo del, by the Poisson-LogNormal mo del, and w e prop ose the Po wer Law mo del as an alternative. In our data, the last mo del seems to b etter describ e the strong heterogeneit y . In order to understand the epidemiological consequences, we use a computational mo del to repro duce a p eculiar w ay of transmission, observed in some cases in nature, where uninfected ticks acquire an infection by feeding on a host where infected ticks are present, without any remark- able epidemiological inv olv ement of the host itself. In particular, w e are in terested in determining the conditions leading to pathogen s pread. W e observ e that the effective transmission of this infection in nature is highly dep enden t on the capabilit y of the im- plemen ted mo del to describ e the tick burden. In addition, we also consider seasonal c hanges in tick aggregation on mice, showing its influence on the spread of the infection. Intro duction Sev eral ecological studies hav e shown that the distribution of tic ks on their hosts is often highly aggregated, with a large num b er of hosts harb ouring few parasites and a small num ber harb ouring a large num b er of them ([14, 15, 34, 2, 25]; other interesting references could be found in [46]). In addition, the distribution of tick developmen t stages is coincident, rather than indep enden t [32]. Sp ecifically , those hosts feeding larv al tick stages were simultaneously feeding the greatest n umber of n ymphs. As a result, ab out 20% of all hosts feed 80% of b oth larv ae and nymphs and the num b er of larv ae feeding alongside n ymphs is t wice as many as it would b e if the distributions were indep enden t [41, 32]. The aggregation of parasites on hosts b ears imp ortan t implications for vector- b orne disease dynamics, since the small fraction of hosts supp orting the bulk of the v ector p opulation is also resp onsible for the ma jorit y of the pathogen transmission [55]. The transmission of tick-borne diseases is characterised by an intricate set of ecolog- ical and epidemiological relationships b et ween pathogen, tic k v ector, v ertebrate hosts and h umans that largely determine their temp oral and spatial dynamics [38]. Tic k- b orne disease dynamics feature several complexities, due to the presence of a n umber of heterogeneities in the system coupled with non-linear phenomena op erating i n the trans- mission pro cesses b et ween ticks, host and pathogen [42]. The transmission of pathogens from one tick to another, a pre-requisite for the establishment of cycles of infection, ma y o ccur via three different path wa ys dep ending on the pathogen (see [36] for a comprehen- siv e review). First, adult female ticks ma y transmit the pathogen to eggs trans-o v arially . Second, tic ks ma y infect a host during their bloo d meal, leading to a systemic infection in the host; tic ks might then acquire the infection b y feeding on an infected host, maintain- ing the infection trans-stadially . Third, ticks may b ecome infected by co-feeding with infected ticks on the same host. Co-feeding transmission is also called non-systemic as it do es not require the host to hav e a systemic infection, since pathogens are transmitted from one tic k to another as they feed in close pro ximity . V ertebrate hosts may v ary in their comp etency to supp ort systemic and co-feeding transmission [24]. Tick-borne pathogens differ also for the mec hanisms which they use to p ersist in nature. F or in- stance, Rickettsia spp., the pathogen agents causing Ro c ky Mountain Spotted F ev er, are 2 main tained by systemic and trans-ov arial transmission in Dermac entor variabili and an- dersoni [1] while it has b een observed that Borr elia Bur gdoferi s.l. spiro c haetes p ersist in nature b y taking adv an tage of all three routes of transmission in I. ricinus , [18, 24]. In the case of the tick-borne encephalitis virus (TBEv), which is an increasing public health concern in Europ e [39, 50, 48], trans-o v arial transmission seems to b e relatively rare and its contribution is generally thought to b e negligible [29]. On the other hand, b oth systemic and non-systemic transmission can take place on reserv oir-comp eten t ro- den t hosts. Ho wev er, due to the v ery short duration of the TBEv infection in ro den ts, [40], the systemic route would only allow infection of a very limited num ber of ticks. Indeed, non-systemic transmission through co-feeding ticks is a more efficient trans- mission route for TBE [40, 41]. Different studies ha ve sho wn that TBEv would not b ecome established in comp eten t hosts, such as ro den ts, without the amplification of the ov erall transmission efficiency provided b y co-feeding transmission (see for instance [40, 45, 20, 19]). The aggregation pattern of tic ks on hosts therefore plays a more im- p ortan t role in the transmission of TBEv than in other tick-borne pathogens, suc h as Borr elia bur gdorferi sensu lato and Anaplasma phago cytophilum, where other efficien t routes of transmission ha ve b een observ ed. Tic k aggregation on hosts and correlation of tick stages facilitate co-feeding transmis- sion and thus significantly increase the basic repro ductiv e n umber, R 0 , of the pathogen, with direct implications for its p ersistence [44, 19]. Using differen t levels of aggregation (from indep enden t to coincident aggregated distribution), Harrison and collab orators [19] show ed that v alues of R 0 increase with progressive levels of aggregation, making it more likely for tic k-b orne pathogens to b ecome established and p ersist. In addition, the authors of the cited w orks evinced that when ticks follow ed a coinciden t aggregated distribution, the increase of R 0 w as greater than in the case of indep enden t aggregated distributions. The degree of aggregation of ticks can b e measured in a n umber of w ays. Since the app earance of influential w orks by Randolph [35] and Shaw et al. ([46] and [47]) the negativ e binomial (NB) distribution has b een extensively used to describ e tick aggrega- tion on hosts (see e.g. [9, 21, 19]). Alternativ ely , other works suggested that different distributions characterised by larger tails than NB (i.e., predicting more ro dents with v ery large tick burden than exp ected with NB), can b e effective in describing tic k aggre- gations. Sp ecifically , a Poisson-LogNormal (PoiLN) mixed mo del has b een successfully used to describ e tick distribution on red grouse chic ks [16], while Bisanzio and collab o- rators [7] show ed the first evidence that the distribution heterogeneity of ticks on hosts seemed to b e b etter describ ed by a p o wer-la w (PL) than a negative binomial distribu- tion. A suitable description of the distribution tail might ha ve important consequences on the dynamics of the pathogen spreading pro cess. Mo delling the spread of vector- b orne diseases through bipartite netw orks [7] show ed that the extreme aggregation of tic ks on hosts has dramatic consequences on the b eha viour of the epidemic threshold. In the current study w e used an extensiv e data set of Ixo des ricinus ticks feeding on mice (a total of 4722 parasitised hosts collected in 9 years) to detect the b est fit for the distribution of tick burden on mice b y testing the p erformance of NB and PoiLN v ersus PL distribution, with particular interest in the shap e of the distribution tail 3 whic h is crucial to suitably describ e the fraction of co-feeding ticks necessary for TBEv transmission. Then, w e used a sto c hastic mo del to simulate the effect of fitting different tic k distributions on the infection dynamics of a tick-borne pathogen. Specifically , we in vestigated the spread of a non-systemically transmitted pathogen (e.g. TBEv) by mo delling the pathogen transmission b et ween susceptible and infectiv e tic ks, considering only co-feeding transmission and distributing tic ks on mice under the hypotheses of NB, P oiLN, and PL distributions. Finally , w e in vestigated the seasonal v ariations in the pattern of tick burden distribution on mice and its implication on TBE-like infection dynamics. Materials and Metho ds Ethics Statement All animal handling pro cedures and ethical issues were approv ed by the Pro vincial Wildlife Managemen t Committee (renewed authorisation n. 595 issued on 04.05.2011) Tick Burden Data Ro den t tick burden data was collected by trapping mice using capture-marking-recapture tec hniques during 2000-2008. The study area was a mixed broadleaf w o o dland [32, 45], lo cated in V alle dei Laghi within the Autonomous Province of T rento, in the north- eastern Italian Alps (grid reference 1652050E 5093750N, altitude 750-800 m a.s.l.). In the year 2000, mice w ere monitored in nine selected areas through placement of 8x8 trapping grids with a 15-m inter-trap in terv al. In 2001 and 2002 the num b er of trapping grids w as reduced to eight, while from 2003 onw ard their num b er was further reduced to four. In summary , the trapping effort consisted of 129 twice-daily trap sessions with at least one capture, resulting in a total num b er of 4722 Ap o demus flavic ol lis captured with at least one tick attached. F or each captured ro den t the n umber and life stage of feeding tic ks was carefully assessed and registered, without remo v al [32, 45]. A total n umber of 55411 ticks w ere counted of which 98 . 64% were larv ae, 1 . 30% were nymphs, and 0 . 04% w ere adults. The num b er of tic ks [nymphs] p er ro den t was b et ween 1 [0] and 111 [15] with a median n umber of ticks p er ro den t equals to 8. Detailed data, on a y early scale, are rep orted in T able 1, while the fraction of n ymphs observed in differen t year and grids is rep orted in T able 2. In Figure 1 the n umber of captured Ap o demus flavic ol lis p er trapping session is sho wn for the whole nine year perio d and for different grids (from A to I ). Data Analysis Tick Burden Distribution Tic ks patterns hav e usually b een describ ed as highly aggregated. Therefore, since the seminal w orks by Crofton [15], Plowrigh t et al. [34], and those b y Anderson and May 4 [2] and [25], the negativ e binomial (NB) probability distribution, q ( k ) =  k + r − 1 k  (1 − p ) r p k , k = 0 , 1 , 2 , . . . (1) has b een considered suitable for describing macroparasite distribution on hosts. Here w e used a maximum-lik elihoo d-estimation (MLE) metho d to estimate the parameters p and r of the probabilit y distribution of the tick burden on the entire dataset obtained b y aggregating capture sessions and grids. In addition, w e considered subsets of the original dataset comp osed b y mice with large num b ers of feeding tic ks to ev aluate the capabilit y of the NB distribution to fit the tail of the parasite distribution. In particular, w e estimated the parameters of the NB distribution on data characterised b y k ≥ k min , where k min represen ts the threshold v alue of ticks p er host ab o ve which the distribution is fitted. T o ev aluate the p erformance of the obtained fits we used Kolmogorov-Smirno v (KS) statistics. The go odness of each fit (GOF) was also ev aluated through a b ootstrap resampling pro cedure, generating 10 3 syn thetic data sets. The obtained p -v alue is de- fined as the relativ e num b er of times that the KS statistic of the fitted distributions on syn thetic data exce eds that measured on real data. Therefore, the larger the p -v alue, the low er our confidence in rejecting the fit. W e considered the conse rv ative v alue 0 . 1 as our threshold v alue, as suggested by Clauset and collab orators [13]. As a first alternative to the NB distribution, we considered a p o wer-la w probability distribution (PL), in its discrete v ersion q ( k ) = Ak − α , k = k min , k min + 1 , k min + 2 , . . . (2) since it may represent a go od candidate to describ e the tail of the distribution [7]. W e recall that A − 1 = P ∞ n =0 1 ( k min + n ) α represen ts the normalising factor of the probability distribution [13]. T o estimate the scaling parameter α of the distribution in suc h a w ay that the PL fits the data for k ≥ k min , we follow ed the algorithm prop osed by Clauset et al. [13]. In short, the fitting pro cedure provides the b est estimate for the parameters k min (called k PL min ) and α by means of MLE and minimisation of KS statistics. F urthermore, b ootstrap techniques w ere used to assess parameter standard deviations (std). W e generated synthetic data and obtained a p -v alue through KS statistics to indicate the go odness of the fit, as for the NB distribution [13]. Another aggregated distribution used for describ e pattern of macroparasites [16] is the P oisson-LogNormal (PoiLN) distribution, q ( k ) = (2 π σ ) − 1 2 k ! Z ∞ 0 λ k − 1 e − λ e − ( log( λ − µ ) 2 ) 2 σ d λ, k = 0 , 1 , 2 , . . . (3) firstly in tro duced b y Bulmer [10] and used in sev eral fields for its capability in describing aggregated data, e.g. [52, 17, 56]. As for the NB distribution, we used a MLE metho d to estimate parameters µ and σ on the entire data set. Uncertain t y on the parameter estimation was assessed b y b o otstrap tec hniques. Moreov er, in line with analysis p er- formed for the NB distribution, w e also explored the capability of the PoiLN distribution 5 to describ e a tick burden larger than a certain threshold k min b y coupling KS statistics and b ootstrap pro cedures. Finally , we compared the PL hypothesis in fitting the tail of the real data distribution with the t w o alternativ es NB and PoiLN by a log-lik eliho o d ratio (LLR) test for differen t v alues of k min . In particular, since the distribution mo dels are non-nested, we used the metho d prop osed b y V uong [51] to understand whether the sign of such test w as statistically significan t or not. Bey ond the estimate of the tic ks-p er-host distribution, we also in vestigated how the tic k burden distributions v ary ov er time and whether a significant difference w as ob- serv ed when different time p erio ds were considered. In particular, w e inv estigated the tic k aggregation patterns during p erio ds characterised by low and high A. flavic ol lis abundance. T o achiev e this goal we smo othed the time series of captured mice with a quadratic p olynomial curv e. The parab ola describing the mice abudance in a sp ecific y ear and grid was normalised b et w een 0 and 1 b efore isolating the time window where this normalised parab ola was higher than a threshold v alue θ ∈ [0 , 1], thus identifying the p eak time of mice abundance, as rep orted in Figure 2. The distribution of ticks feeding on mice has b een ev aluated and compared considering in- and out-of-p eak time p eriods for differen t v alues of θ . W e calculated the KS statistic b et ween the in- (high abundance) and out-of- (low abundance) p eak time distributions of tic k burden, and we then compared the v alue observed in real data to a b o otstrapped data set in order to es- tablish whether this measure was statistically significant. F or this purp ose we generated 10 5 syn thetic in-and out-of- p eak samples having the same size as the observ ed ones. As a test of soundness, w e then calculated the fraction of the KS statistic that is larger in syn thetic data than on real data. La rval and Nymphal Aggregations P atterns on Mice A necessary condition for an effectiv e non-systemic transmission of a pathogen is the co- incidence of the larv al and nymphal aggregation distributions on hosts [41]. Therefore, our first step was to examine the asso ciation b et ween the n umber of larv ae and that of n ymphs on each host using Sp earman’s rank-correlation co efficien t [57]. In particular, w e preferred a non-parametric metho d rather than the more commonly used Pearson’s correlation co efficien t since tick distributions are aggregated (i.e. deviate from normal distribution) and we were more in terested in any monotonic relations of our v ariables than in the linear relation depicted by the Pearson’s co efficien t. More in detail, a p os- itiv e [negative] Sp earman’s co efficient would indicate that an increase in the num b er of nymphs p er mouse is asso ciated with an increase [decrease] of the num b er of larv ae p er mouse. Therefore, a p ositiv e Sp earman’s correlation co efficien t could b e interpreted as an indicator of the coincidence of the distributions, a zero co efficien t could suggest the indep endence of the tw o distributions, and a negative co efficien t, an uncommon re- sult, w ould b e an indicator of ha ving tw o unimo dal distributions with tw o async hronous p eaks. Moreov er, to ev aluate the significance of Spearman’s coefficient (i.e. the probabil- it y that the same co efficien t could be obtained b y c hance) w e implemen ted a p erm utation test. In particular, w e compared the ev aluations on synthetic datasets with a resh uffled 6 n umber of nymphs and on the original data and counted the num ber of times that the absolute v alue of Sp earman’s co efficien t w as larger than for the original data. The lo wer the sum, the higher our confidence in in terpreting the asso ciation as significan t. T o further ev aluate the coincidence of tic k stage distributions and the consequences on the non-systemic transmission of a pathogen, following Randolph et al. [41], w e ev aluated the mean num b er of larv ae cofeeding with a n ymph on a host. In fact, the larger the mean, the larger the num b er of larv ae that can p otentially b e infected via non-systemic transmission. After obtaining this empirical datum, w e calculated the mean v alue for 10 3 syn thetic datasets where the num b er of nymphs w as reshuffled, simulating indep enden t distributions in order to hav e a more robust interpretation. After comparison of empirical and synthetic datasets, a significan tly larger empirical mean n umber of larv ae p er n ymphs giv es evidences of coincident distributions, [41]. Simulations of Tick-Borne Disease Spreading via Non-Systemic T ransmission In order to explore the impact of different parasite aggregation distributions on the spread of a TBEv-lik e pathogen where the main transmission route is through co-feeding, w e p erformed extensive numerical simulations informed by the data ab out tick aggre- gation on mice. In this setting, tic k larv ae were not infective (transov aric transmission has b een indicated as negligible [22]), adults only rarely feed on mice (on our data set adults ticks are ab out 0 . 05% of the total num b er of tic ks feeding on mice), and the only transmission link that w e considered was the co-feeding b et w een infectiv e nymphs and larv ae. Therefore, the only actors in our mo del w ere n ymphs and larv ae feeding on hosts. Moreov er, Ros` a and collab orators suggested in a recen t w ork devoted to the same geographical area [45] that the larv ae that feed in one y ear generally quest and feed as nymphs in the following year. Therefore, by adapting the Susceptible-Infected- Susceptible (SIS) mo del [3] to our purp ose we assumed that nymphs are categorised as infectiv e or not, that feeding larv ae are susceptible and that some of them could even- tually b e infected by co-feeding with infective n ymphs b efore moulting (thus b ecoming infectiv e n ymphs at time t + 1). A t each iteration t , with t b eing a discrete n umber b et w een t 0 and t max and ∆ t = 1 year, w e assigned a num b er of ticks to each of the N h mice by drawing a sample from the considered distribution q . Then, on each mouse w e said that of k tic ks feeding on it, k f were n ymphs and the other larv ae (with 0 < f < 1). These nymphs were larv ae in the previous year and were p ossibly infected. Then, defin- ing as π L ( t − 1) the prev alence among larv ae after feeding at time t − 1, w e assumed that the prev alence at time t among n ymphs was π N ( t ) = π L ( t − 1). Th us, the num b er of infectiv e nymphs on a mouse that at time t was parasitised by k ticks w as k f π N ( t ). Then, on eac h of the N h mice the co-feeding transmission b et ween larv ae and infective n ymphs could o ccur with probability β and w e up dated π L ( t ) accordingly to the fraction of larv ae infected (i.e. the fraction of infective n ymphs at next time step). The follo wing meta-co de summarises the epidemiological dynamic 1. for t b et ween t 0 and t max : 7 a) for each mouse i , with i b et ween 1 and N h • k ( i ) is the num b er of ticks it feeds, b eing k ( i ) a n umber drawn from the probabilit y distribution q • of the k ( i ) tic ks, f k ( i ) are n ymphs and the remaining larv ae • of the f k ( i ) nymphs, a fraction π L ( t − 1) f k ( i ) are infective, the others are susceptible • non-systemic transmission b et ween infective nymphs and larv ae on the same host o ccurs with probability β b) π L ( t ) is up dated as the fraction of larv ae infected c) if π L ( t ) is equal to zero w e stop the lo op It is worth stressing that in the previous meta-co de we did not consider ticks reco vering from the infection, since we assumed that a feeding infective nymph at time t will exit the infectious dynamics b y moulting to the adult stage or dying. W e also mo dified the previous dynamics to deal with different distributions in tic k aggregation as a function of seasonality . At eac h year t , we classified mice as observed during the mice p eak activit y (= γ N h mice, with 0 < γ < 1) and observed out of the p eak (= (1 − γ ) N h ). Therefore, w e assigned the num ber of ticks feeding on mice according to the respectively aggregated distributions q IN and q OUT . Moreov er, since the larv ae obtaining a blo od meal at year t will b e nymphs at year t + 1 without an y other in volv emen t in the epidemic spreading at y ear t , [45], these mo difications to the meta- co de are sufficien t to suitably describ e the seasonal v ariation in the epidemic pro cess. More explicitly , the epidemic dynamic in the presence of seasonality in tick aggregation ma y b e describ ed by the following meta-co de: 1. for t b et ween t 0 and t max : a) a fraction γ of the N h mice are lab elled as observed during mice p eak activity (the remaining = (1 − γ ) N h as observ ed out of the p eak windo w) b) for each mouse i , with i b et ween 1 and N h • k ( i ) is the num b er of ticks it feeds, b eing k ( i ) a n umber drawn from the probabilit y distribution q IN , if the mouse was lab elled as observed during the mice p eak activity , or q OUT , if not • of the k ( i ) tic ks, f k ( i ) are n ymphs, the remaining larv ae • of the f k ( i ) nymphs, a fraction π L ( t − 1) f k ( i ) are infectiv e, the other susceptible • non-systemic transmission b et ween infective nymphs and larv ae on the same host o ccurs with probability β c) π L ( t ) is up dated as the fraction of larv ae infected d) if π L ( t ) is equal to zero w e stop the lo op 8 Results Ticks Burdens The probabilit y distribution of tick burden on mice was skew ed and sho wed a heavy tail. The b est fit of the NB distribution was obtained on the largest a v ailable subsets of data, i.e. with k min = k NB min = 1, see left panel of Figure 3. In this setting, the MLE metho d estimated r = 1 . 30 (95% confidence in terv als (CI) = 1 . 25 , 1 . 35) and p = 0 . 10 (95%CI=0 . 09 , 0 . 10). How ever, the GOF of the NB distribution was very low ( p < 10 − 3 ) for any v alue of k min , see central panel of Figure 3, thus giving evidence for rejecting the h yp othesis of the NB functional form. Similarly , the b est fit of P oiLN distribution w as ac hieved on the largest subsets of data, ( k min = k PoiLN min = 1, see left panel of Figure 3). In this case the estimated parameters w ere µ = 1 . 96 (95%CI=1 . 92 , 1 . 99) and σ = 0 . 99 (95%CI=0 . 96 , 1 . 02). The GOF of the P oiLN, central panel of Figure 3, suggested that P oiLN was acceptable only for k min > 38. How ev er, for k min > 38, the KS statistic displa yed v alues that were to o large to consider the P oiLN distribution appropriate for describing real data. On the other hand, by fitting the tail of the distribution to a PL distribution, we found that the b est fit w as obtained for k min = k PL min = 38 (with a standard deviation of 5 . 83), see left panel of Figure 3. This k min v alue is matched with an estimated scaling parameter α = 4 . 27 (with standard deviation = 0 . 41). The GOF test ( p -v alue larger than 0 . 1) suggested that the optimum PL fit on the tail of the distribution should not b e ruled out, and that the result holds for every PL fit with k min > 35 see center panel of Figure 3. Finally , the LLR test highlighted that the PL fitting is to b e preferred ( p < 10 − 3 ) to the NB in describing the tail of the distribution for a large range of low er b ounds, k min ∈ [8 , 44], see right panel of Figure 3. Similarly , the PL is to b e preferred to the P oiLN for k min ∈ [5 , 54]. Moreov er, it is worth to stress that for v alues ab o ve 44 (55) the sign of the LLR test still indicates the PL fit as the preferred one compared to the NB (and P oiLN), although the indication loses statistical significance due to the scarcit y of av ailable data. In Figure 4 w e sho w the complementary cumulativ e probability distribution of the best fits resulting from k NB min = k PoiLN min = 1 for NB and PoiLN distributions and k PL min = 38 for PL distribution against field data of the n umber of ticks p er mouse. F rom this plot w e noticed that ab o ve a certain n umber of ticks p er mouse NB [P oiLN] under-estimates [o ver-estimates] the tail of the distribution (indeed b oth fits were statistically ev aluated as very p o or). A t the same time, in agreement with statistical results summarised in Figure 3, we noticed that the PL fit in Figure 4 more appropriately describ es the righ t tail of the data distribution. The num b er of mice captured in different y ears and grids show ed strong seasonal patterns as rep orted in Figure 1. F or each grid and each year we defined tw o separate p eriods dep ending on the mice abundance as defined in section “Data Analysis” and sk etched in Figure 2. Imposing a threshold θ , for eac h year and grid we identified a time window of high mice abundance. With θ = 0 . 5 we found significant evidence that the distribution of ticks on mice within the abundance p eak was differen t from that 9 observ ed outside. Indeed, the fraction of the KS measures calculated on the synthetic samples lo wer than the real-data KS statistic was almost 98%, th us indicating very low confidence in obtaining the same measurement by c hance. The same statistical evidence w as also obtained by using different time window thresholds (such as, θ = 0 . 4 and 0 . 6). On the data sets classified as inside (IN) and outside (OUT) the time windo w of mice abundance p eak, w e fitted for different time-window lengths ( θ = 0 . 4 , 0 . 5 , 0 . 6) the parameters r and p for NB distribution (Figure 5, left panels) and α and k min for PL distribution (Figure 5, right panels). W e observ ed a larger PL scaling parameter α inside the mice abundance p eak than outside (tw o-sample t -test output: for θ = 0 . 5 t -statistic= − 74 . 95, df= 1931, p < 10 − 3 ) indicating a larger heterogeneity in tick burden outside the abundance p eak time. Moreov er the GOF test indicated a rejection of the NB fit in b oth sets (IN and OUT) with θ = 0 . 5. On the other hand, the GOF test with θ = 0 . 5 show ed that the PL mo del cannot b e ruled out in b oth sets ( p -v alue > 0 . 1) and the LLR test indicated that the PL fitting outp erforms the NB mo del ( p -v alue < 0 . 05) in the estimates b oth inside and outside the p eak time window. The distribution of larv ae and nymphs on mice are coinciden t rather than indep enden t, and indeed the same 20% most infested hosts feed b oth 55% of the nymphs and 54% of the larv ae. Moreov er, Sp earman’s correlation co efficien t measured on the num ber of larv ae and nymphs on mice was p ositiv e (0 . 24) and the probabilit y that this co efficien t w as detected by chance w as very low (the empirical v alue was the largest if compared to those ev aluated in 10 3 resh uffled samples). In addition, the mean num b er of larv ae co-feeding with a n ymph is ab out 23 which is almost double the v alue that w ould b e seen if the distributions w ere indep enden t (mean equal to 12). Non-systemic disease spreading simulations T o start, we simulated the non-systemic disease spreading of a TBE-lik e pathogen with a fraction f of nymphs among ticks equals to 2%, close to the one observed in our real data (cfr. T able 2), 5%, and 10%, as in literature [37, 41]. W e consider the empirical distribution observ ed on the entire data set. W e fixed the num b er of hosts to N h = 10 4 whic h, together with the considered distribution, resulted in a n umber of vectors pairs equal to N V ∼ 10 5 . In our simulations, we explored the effects of β , the infection probabilit y , on the observed prev alence at the final time step, π L ( t max ), with t max = 1000. (W e observ ed that t max = 1000 was larger enough to allo w the prev alence to conv erge to ward an endemic pseudo-equilibrium or the disease-free equilibrium). F or each β we allo wed 200 simulations to run starting from an initial prev alence of π N ( t 0 ) = 1%. In Figure 6 w e plotted the prev alences (median v alue, interquartile in terv als and the 95%CI) observ ed at equilibrium as a function of the transmission probabilities, β . Results sho wed that the larger the fraction of nymphs among ticks feeding on mice, the larger the probabilit y of pathogen inv asion and the infection prev alence. Then, w e explored the effects of differen t tic k burden distributions on the spread of infection. T o this end we considered four distributions: PL, NB, PoiLN, and the empirical distribution on the entire data set (aggregated on capture sessions and grids). F or synthetic distributions w e considered the actual observed distribution b elo w the 10 estimated k min , while we used the b est fit of synthetic distributions to describ e v alues greater than k min . Again, we fixed the n umber of hosts to N h = 10 4 . It is w orth stressing that in the syn thetic samples generated from these distributions we observed some features similar to those observ ed in real-data. F or instance, the num ber of nymphs w as p ositiv ely asso ciated with that of larv ae and more particularly a nymph co-fed with a mean num b er of larv ae similar to that observ ed in realit y (for PL the mean num b er w as 23, for NB 20, and for PoiLN 27). Results, plotted in Figure 7 for f = 2% and in S2-text for f = 5% and f = 10%, corrob orated the hypothesis that the transmission probability needed for the pathogen to b ecome endemic is driven by the shap e of the tail of the distributions. In particular, w e noticed that for the PoiLN distribution (the one with larger fitted tail) the epidemic threshold is the low est, while for the NB distribution (the one with smaller fitted tail) the infection probabilit y needed for in v asion is the highest. Not surprisingly , the PL, which has the b est p erformances in fitting the tail of the empirical distribution, is the one for whic h the prev alences at equilibria b etter resem ble those observed in simulations using the empirical distribution. W e also p erformed some sensitivity analysis on parameter distributions, further highlighting that the larger the tail of the distribution, the lo wer the epidemic threshold (see SM.1). In addition, sensitivity analysis on the fraction of n ymphs ( f ) show ed that f do es not qualitatively influence the epidemic b ehaviour (see SM.2). F urthermore, we inv estigated the effect of differences in the distribution of the tic k burden as a function of the abundance of mice on the spreading of a non-systemic infectious disease. T o this end, w e fixed γ = 0 . 89, as measured in the dataset, and as q IN w e considered a PL with exp onen t α IN = 4 . 39 as estimated with θ = 0 . 5. In a similar wa y , we assumed as q OUT a PL distribution with exp onen t α IN = 3 . 48. F or b oth q IN and q OUT w e further set k min = 5 , 10 , 15. Results are summarised in Figure 8, from whic h it could b e inferred that the epidemic outcome was strongly influenced b y the different distributions of feeding tic ks according to mice abundance. W e consistently observ ed that the transmission probability needed for the pathogen to effectively spread w as smaller when the time windows identified by mice abundance are considered. Discussion Tic k aggregation on hosts is the result of several complex interactions of biotic and abi- otic factors, such as host exp osure and susceptibility to tic ks, tic ks’ phenology and host b eha viour, environmen tal factors, av ailability of resources, and others [53, 9]. Histor- ically , the NB distribution has b een preferred to the P oisson distribution to describ e parasite heterogeneity across hosts b ecause it suitably repro duces ov erdispersed obser- v ations. It has also b een widely used in empirical [35, 46, 47, 21] and theoretical studies [44, 43, 19]. How ev er, fat taile d distributions other than the NB one can also adequately repro duce tic k aggregation, as shown by Elston et al. [16] and Bisanzio and collab orators [7]. Through the use of an extensive data set of feeding Ixo des ricinus ticks on mice, 11 w e show ed that a PL distribution is b etter able to describ e the right tail of the tick distribution on hosts than a NB or a P oiLN distribution (see Figure 3 and 4). This finding may hav e relev an t epidemiological consequences, since it is well do cumen ted that the heterogeneity of contact distributions among individuals has large impacts on pathogen spread and p ersistence [6, 4, 23, 30, 27, 26, 5]. In fact, it has b een demonstrated [31] that the minimum transmission probability for a pathogen to spread on a net work, the so-called epidemic threshold, is driv en b y the first and the second moment of this distribution. In particular, P astor-Satorras et al. [31] demonstrated that the larger the heterogeneit y , the low er the epidemic threshold for the pathogen to spread, with an interesting b eha viour in infinite size net work sho wing a zero epidemic threshold [30]. Th us, the epidemiological inferences on the spread of a pathogen are highly influenced by the c haracterisation of the connectivit y distribution and in particular b y the distribution tail (i.e. the heterogeneit y). Our results corrob orate those findings and generalise them in a different framew ork and for more complex transmission routes, i.e. a vector-host net work for non-systemically transmitted diseases. In particular, w e found that the tail of the distribution of the num b er of tic ks p er ro den t highly influences pathogen spreading (see Figure 7 and S1-text). F urthermore, it is w orth remarking that although the tail of the distribution as defined here represents ab out 5% of the entire data set, our simulation findings suggest that this small part of the distribution is crucial for pathogen inv asion. W e also confirm that the probability of pathogen in v asion and the infection prev alence are strongly influenced by the fraction f of nymphs on the total feeding tic ks on mice (Figure 6 and S2-text). The co-o ccurrence of larv ae and n ymphs on comp eten t hosts is in fact essential for the horizon tal transmission of non-systemic transmitted tick-borne pathogens, suc h as TBE, and it has b een do cumented, b oth empirically and theoretically , that it could b e a key factor in creating TBE hotsp ots, [11, 8]. Our conclusions confirm previous findings sho wing that the distribution of ticks on ro den ts ma y significantly affect the spread of infections [9, 7, 12], esp ecially for non- viraemic transmitted diseases such as TBE [44, 32, 19]. Under the h yp othesis of a NB distribution of ticks across hosts, b oth Ros` a et al. [44] and Harrison and collab orators [19] sho w ed that highly coincident and aggregated distributions fa vour the establishmen t of TBEv. Ho wev er, highly heterogeneous degree distributions do not necessary imply a higher spread of disease. Indeed, Piccardi et al. [33] show ed that scale-free netw orks can b e muc h less efficien t than homogeneous netw orks in fav ouring the disease spread in the case of a nonlinear force of infection. The correct description of tic k aggregation on hosts could dramatically affect disease con trol strategies: for instance, P erkins [32] emphasised that an optimised con trol effort targeted on highly parasitised mice, also identified as sexually mature males of high b ody mass, could significan tly low er the transmission p otential. On the other hand, Brunner and colleagues [9] observ ed that the identification of individuals whic h fed a disprop ortionate num b er of ticks (and that can therefore act as sup erspreaders) can b e c hallenging, since simple cov ariates such as sex, age or mass do not en tirely explain the differences in parasite burden. In order to fully understand the differen t tick attachmen t b eha viours on hosts, we iden tified different time windows related to ro den t seasonal dynamics. Using this ap- 12 proac h w e found that the distribution of ticks on mice ma y v ary across the season, with higher aggregation heterogeneity in p erio ds of lo w ro den t abundance and low er aggregation heterogeneity during the p eak of host abundance (see Figure 5). W e also sho wed that seasonal aggregation patterns, characterised by larger tails in time perio ds of low host abundance, enhance the spread of non-viraemic transmitted diseases (see Figure 8). Shaw and collab orators [47] observed significant v ariations in the degree of aggregation b et w een host subsets – stratified by sex, age, space or time of sampling – in several host-parasite systems. In agreemen t with our results (low er aggregation in p eriod of high mice abundances as shown b y estimated exp onen ts of PL), they found that aggregation in cop epo d ( L ep e ophtheirus p e ctor alis ) infesting plaice ( Pleur one ctus platessa ) decreases during summer months. They mainly ascrib ed the observed v ari- ation to significant differences in mean parasite burden among months. On the other hand, we did not find significant differences in tick burden inside and outside the win- do w of high ro den t abundance. Sp ecifically , in the case of θ = 0 . 4 , 0 . 5 , 0 . 6, the a verage n umber of tic ks p er host were 11 . 96 , 11 . 96 , 11 . 84 inside the window of high ro dent abun- dance and 12 . 34 , 12 . 23 , 12 . 78 outside and the differences b et ween inside and outside are not statistically significant (p erm utation tests, p > 0 . 05). Ho wev er, the second moment of the num b er of ticks p er host drastically c hanged b et w een high and lo w abundance p eriods, driving the difference in the aggregation distributions observed in the t wo time windo ws. Seasonal v ariations in resource av ailability and host abundance can hav e a significan t effect on the space used by mice. Males and females tend to resp ond to these c hanges in different wa ys, since space use for females is driven largely by fo o d av ail- abilit y , whereas the distribution of males is related primarily to mating opp ortunities. Y ello w-neck ed mouse ( A. flavic ol lis ) females exhibited reduced spatial exclusivity and larger home ranges during low er fo od av ailability while males v aried their spatial dis- tribution accordingly by also expanding their home ranges [49]. An inv erse relationship b et w een p opulation densit y and home range sizes has also b een observed in woo d mice ( Ap o demus sylvaticus ) [54]. Consequen tly , in p eriods of low rodent abundance more mo- bile ro den ts, esp ecially males, are more likely to hit a patch of larv al tic ks. As a result, these individuals w ould harb our a large amoun t of ticks and increase the aggregation of tick distribution among the ro den t p opulation. On the other hand, tic k density is usually low er in p erio ds of low ro den t abundance, and the av erage tick burden would decrease for the rest of the p opulation, esp ecially females, balancing the ov erall tic k bur- den. On the contrary , during times of high abundance mice mov e less and tic ks would b e distributed more evenly among the ro den t p opulation resulting in the observ ation of a lo wer aggregation in tick distribution during the p eak of ro den t abundance. Our primary goal was to help understand the role of tic k aggregation across mice on the spread of non-viraemic transmitted diseases through a simple and general transmission mo del. Other works – such as [44, 28, 43, 8] – describ ed in v ery fine detail the trans- mission of v ector-b orne diseases, in tro ducing differen t transmission routes, tick stages and alternative hosts in the epidemic mo del. F or instance, Norman and colleagues [28] demonstrated through an epidemiological mo del that non-viraemic transmission could ha ve non-negligible effects on the p ersistence of a disease lik e the Louping ill. Here, considering the non-systemic transmission only , we explored the effect of using different 13 theoretical functional forms to describ e the tick burden on hosts. By estimating param- eters of the burden distributions on a very detailed data set, we defined a simple and transparen t transmission mo del that explicitly takes in to account the real contact pattern of vectors and hosts in the description of a non-system atically transmitted v ector-b orne disease. In this w ay w e were able to emphasise that, while the NB and PoiLN mo dels can sufficiently fit the whole real distribution, the PL mo del represents a b etter fit for the distribution tail. F urthermore, the vector p ersp ectiv e approach used in our mo del giv es b etter insights into the dynamics of non-systemic transmitted pathogens resp ect to host p ersp ectiv e mo dels that w ere more commonly and widely used in this context [44, 43, 28, 8]. In addition, epidemiological simulations parameterised by the fitted tick burden distributions highligh ted the epidemiological consequences of describing tic k ag- gregation on hosts trough distributions with different tails, sho wing that the shap e of the tail distribution has a non-negligible influence on pathogen p ersistence. F uture works will b e devoted to extend the present findings to more complex transmission dynamics (e.g. including viraemic or transov aric transmission), in order to assess the effect of a PL deca y of the distribution for a wider range of vector-borne diseases. Ackno wledgements Authors wish to thank all the field assistants from the Department of Bio diversit y and Molecular Ecology at F ondazione Edm und Mach. Authors wish to thank Bryan N. Iotti for the painstaking pro ofreading of the whole manuscript. Financial Disclosure This study was partially funded by Europ ean Union (EU) Grant FP7-261504 EDENext and is catalogued by the EDENext Steering Committee as EDENext265. The conten ts of this publication are th e sole resp onsibilit y of the authors and do not necessarily reflect the views of the Europ ean Commission. LF ac knowledges supp ort from the Lagrange Pro ject, CR T and ISI F oundation. References [1] J.R. Amsden, S. W armack, and P .O. Gubbins. Tick-borne bacterial, rick ettsial, spiro c hetal, and protozoal infectious diseases in the united states: A comprehensive review. Pharmac other apy , 25(2):191–210, 2005. doi:10.1592/phco.25.2.191.56948. [2] R.M. Anderson and R.M. Ma y . Regulation and stability of host-parasite population in teractions: I. Regulatory pro cesses. Journal of Animal Ec olo gy , 47(1):219–247, 1978. [3] R.M. Anderson and R.M. Ma y . Infe ctious Dise ases of Humans: Dynamics and Contr ol . OUP Oxford, 1992. 14 [4] F. Ball, D. Mollison, and G. Scalia-T omba. Epidemics with tw o lev els of mixing. A nnals of Applie d Pr ob ability , 7(1):46–89, 1997. [5] Alain Barrat, Marc Barth´ elem y , and Alessandro V espignani. Dynamic al Pr o c esses on Complex Networks . Cambridge Universit y Press, 2008. [6] Niels Beck er and Ian Marschner. The effect of heterogeneity on the spread of dis- ease. In Jean-Pierre Gabriel, Claude Lefvre, and Philipp e Picard, e ditors, Sto chas- tic Pr o c esses in Epidemic The ory , v olume 86 of L e ctur e Notes in Biomathemat- ics , pages 90–103. Springer Berlin Heidelb erg, 1990. ISBN 978-3-540-52571-4. doi:10.1007/978-3-662-10067-7-9. [7] D. Bisanzio, L. Bertolotti, L. T omassone, G. Amore, C. Ragagli, A. Mannelli, M. Gi- acobini, and P . Pro vero. Mo deling the spread of v ector-b orne diseases on bipartite net works. PL oS ONE , 5(11), 2010. doi:10.1371/journal.p one.0013796. [8] L. Bolzoni, R. Ros` a, F. Cagnacci, and A. Rizzoli. Effect of deer density on tick infestation of ro dents and the hazard of tick-borne encephalitis. I I: Population and infection mo dels. International Journal for Par asitolo gy , 42(4):373–381, 2012. doi:10.1016/j.ijpara.2012.02.006. [9] J.L. Brunner and R.S. Ostfeld. Multiple causes of v ariable tick burdens on small- mammal hosts. Ec olo gy , 89(8):2259–2272, 2008. doi:10.1890/07-0665.1. [10] M. G. Bulmer. On fitting the p oisson lognormal distribution to sp ecies-abundance data. Biometrics , 30(1):101–110, 1974. [11] F. Cagnacci, L. Bolzoni, R. Ros` a, G. Carpi, H.C. Hauffe, M. V alen t, V. T agliapi- etra, M. Kazimirov a, J. Ko ci, M. Stank o, M. Luk an, H. Hen ttonen, and A. Rizzoli. Effects of deer densit y on tick infestation of ro den ts and the hazard of tick-borne encephalitis. I: Empirical assessmen t. International Journal for Par asitolo gy , 42(4): 365–372, 2012. doi:10.1016/j.ijpara.2012.02.012. [12] J.M. Calabrese, J.L. Brunner, and R.S. Ostfeld. P artitioning the aggre- gation of parasites on hosts into intrinsic and extrinsic comp onen ts via an extended p oisson-gamma mixture mo del. PL oS ONE , 6(12):e29215, 2011. doi:10.1371/journal.p one.0029215. [13] A. Clauset, C.R. Shalizi, and M.E.J. Newman. Po w er-la w distributions in empirical data. SIAM R eview , 51(4):661–703, 2009. doi:10.1137/070710111. [14] H.D. Crofton. Nemato de parasite p opulations in sheep on lo wland farms. VI. Sheep b eha viour and nemato de infections. Par asitolo gy , 48(3-4):251–260, 1958. [15] H.D. Crofton. A mo del of hostparasite relationships. Par asitolo gy , 63:343–364, 1971. doi:10.1017/S0031182000079890. 15 [16] D.A. Elston, R. Moss, T. Boulinier, C. Arrowsmith, and X. Lambin. Analysis of aggregation, a work ed example: Numbers of tic ks on red grouse chic ks. Par asitolo gy , 122(5):563–569, 2001. doi:10.1017/S0031182001007740. [17] U. Gonzales-Barron and F. Butler. A comparison b et w een the discrete p oisson- gamma and p oisson-lognormal distributions to c haracterise microbial counts in fo ods. F o o d Contr ol , 22(8):1279–1286, 2011. doi:10.1016/j.fo odcont.2011.01.029. [18] J. Gra y , O. Kahl, R.S. Lane, and G. Stanel, editors. Lyme Bor el liosis - Biolo gy, Epidemiolo gy and Contr ol . CABI publishing, 2002. [19] A. Harrison and N.C. Bennett. The imp ortance of the aggregation of ticks on small mammal hosts for the establishment and p ersistence of tic k-b orne pathogens: An in vestigation using the R 0 mo del. Par asitolo gy , 139(12):1605–1613, 2012. doi:10.1017/S0031182012000893. [20] N.A. Hartemink, S.E. Randolph, S.A. Davis, and J.A.P . Heesterb eek. The basic repro duction n umber for complex disease systems: Defining R 0 for tick-borne in- fections. Americ an Natur alist , 171(6):743–754, 2008. doi:10.1086/587530. [21] C. Kiffner, T. V or, P . Hagedorn, M. Niedrig, and F. Rhe. F actors affecting patterns of tick parasitism on forest ro den ts in tick-borne encephalitis risk areas, german y . Par asitolo gy R ese ar ch , 108(2):323–335, 2011. doi:10.1007/s00436-010-2065-x. [22] L. Lindquist and O. V apalahti. Tick-borne encephalitis. The L anc et , 371(9627): 1861–1871, 2008. doi:10.1016/S0140-6736(08)60800-4. [23] A.L. Lloyd and R.M. May . How viruses spread among computers and p eople. Sci- enc e , 292(5520):1316–1317, 2001. [24] A. Mannelli, L. Bertolotti, L. Ge rn, and J. Gra y . Ecology of Borr elia bur gdorferi sensu lato in Europ e: T ransmission dynamics in m ulti-host systems, influence of molecular pro cesses and effects of climate change. FEMS Micr obiolo gy R eviews , 36 (4):837–861, 2012. doi:10.1111/j.1574-6976.2011.00312.x. [25] Rob ert M. May and Roy M. Anderson. Regulation and stability of host-parasite p opulation interactions: I I. Destabilizing pro cesses. Journal of Animal Ec olo gy , 47 (1):249–267, 1978. [26] Rob ert M. Ma y and Alun L. Lloyd. Infection dynamics on scale-free net works. Physic al R eview E , 64:066112, Nov 2001. doi:10.1103/PhysRevE.64.066112. [27] M. E. J. Newman. Spread of epidemic disease on netw orks. Physic al R eview E , 66: 016128, Jul 2002. doi:10.1103/PhysRevE.66.016128. [28] R. Norman, D. Ross, M.K. Laurenson, and P .J. Hudson. The role of non-viraemic transmission on the p ersistence and dynamics of a tic k b orne virus - louping ill in red grouse ( L agopus lagopus sc oticus ) and moun tain hares ( L epus timidus ). Journal of Mathematic al Biolo gy , 48(2):119–134, 2004. doi:10.1007/s00285-002-0183-5. 16 [29] P .A. Nuttall and M. Labuda. Dynamics of infection in tick vectors and at the tic k- host in terface. A dvanc es in Virus R ese ar ch , 60:233–272, 2003. doi:10.1016/S0065- 3527(03)60007-2. [30] R. P astor-Satorras and A. V espignani. Epidemic spreading in scale- free net works. Physic al R eview L etters , 86(14):3200–3203, 2001. doi:10.1103/Ph ysRevLett.86.3200. [31] Romualdo Pastor-Satorras and Alessandro V espignani. Epidemic dynamics in finite size scale-free netw orks. Physic al R eview E , 65:035108, Mar 2002. doi:10.1103/Ph ysRevE.65.035108. [32] S.E. Perkins, I.M. Cattadori, V. T agliapietra, A.P . Rizzoli, and P .J. Hudson. Em- pirical evidence for key hosts in p ersistence of a tick-borne disease. International Journal for Par asitolo gy , 33(9):909–917, 2003. doi:10.1016/S0020-7519(03)00128-0. [33] C. Piccardi and R. Casagrandi. Inefficient epidemic spreading in scale-free netw orks. Physic al R eview E , 77(2):026113, 2008. doi:10.1103/PhysRevE.77.026113. [34] R.C. Plowrigh t and J.E. Paloheimo. A theoretical study of p opulation dynamics in the sheep tic k. The or etic al Population Biolo gy , 12(3):286–297, 1977. [35] Sarah E. Randolph. Patterns of distribution of the tic k Ixo des triangulic eps birula on its hosts. Journal of Animal Ec olo gy , 44(2):451–474, 1975. [36] S.E. Randolph. Tic ks are not insects: Consequences of contrasting v ector biology for transmission p oten tial. Par asitolo gy T o day , 14(5):186–192, 1998. doi:10.1016/S0169-4758(98)01224-1. [37] S.E. Randolph and N.G. Craine. General framework for comparative quan titative studies on transmission of tic k-b orne diseases using lyme b orreliosis in Europ e as an example. Journal of Me dic al Entomolo gy , 32(6):765–777, 1995. [38] S.E. Randolph and EDEN-TBD sub-pro ject team. Human activities predominate in determining c hanging incidence of tick-borne encephalitis in Europ e. Eur o surveil- lanc e : bul letin eur op en sur les maladies tr ansmissibles = Eur op e an c ommunic able dise ase bul letin , 15(27):24–31, 2010. [39] S.E Randolph and D. Sumilo. Tic k-b orne encephalitis in Europ e: dynamics of c hanging risk. In W. T akken and B. Knols, editors, Emer ging Pests and V e ctor- b orne Dise ases in Eur op e , pages 187–206. W ageningen Academic, 2007. [40] S.E. Randolph, L. Gern, and P .A. Nuttall. Co-feeding tic ks: Epidemiological signif- icance for tick-borne pathogen transmission. Par asitolo gy T o day , 12(12):472–479, 1996. doi:10.1016/S0169-4758(96)10072-7. 17 [41] S.E. Randolph, D. Miklisov, J. Lysy , D.J. Rogers, and M. Labuda. Inci- dence from coincidence: P atterns of tick infestations on ro den ts facilitate trans- mission of tic k-b orne encephalitis virus. Par asitolo gy , 118(2):177–186, 1999. doi:10.1017/S0031182098003643. [42] S.E. Randolph, C. Chemini, C. F urlanello, C. Genchi, R. A. Hails, P .J. Hudson, L.D. Jones, G. Medley , R. Norman, A. Rizzoli, G. Smith, and M.E.J. W o olhouse. The ecology of tick-borne infections in wildlife reservoirs. In P .J. Hudson, A. Rizzoli, B.T. Grenfell, H. Hesterb eek, and A.P . Dobson, editors, The Ec olo gy of Wild life Dise ases , pages 119–138. Oxford Universit y Press, 2002. [43] R. Ros` a and A. Pugliese. Effects of tick p opulation dynamics and host densities on the p ersistence of tick-borne infections. Mathematic al Bioscienc es , 208(1):216–240, 2007. doi:10.1016/j.mbs.2006.10.002. [44] R. Ros` a, A. Pugliese, R. Norman, and P .J. Hudson. Thresholds for disease p er- sistence in mo dels for tic k-b orne infections including non-viraemic transmission, extended feeding and tick aggregation. Journal of The or etic al Biolo gy , 224(3):359– 376, 2003. doi:10.1016/S0022-5193(03)00173-5. [45] R. Ros` a, A. Pugliese, M. Ghosh, S.E. P erkins, and A. Rizzoli. T emp oral v ariation of Ixo des ricinus intensit y on the ro den t host Ap o demus flavic ol lis in relation to lo cal climate and host dynamics. V e ctor-Borne and Zo onotic Dise ases , 7(3):285–295, 2007. doi:10.1089/vbz.2006.0607. [46] D.J. Shaw and A.P . Dobson. Patterns of macroparasite abundance and aggregation in wildlife p opulations: A quantitativ e review. Par asitolo gy , 111(SUPPL.):S111– S133, 1995. [47] D.J. Shaw, B.T. Grenfell, and A.P . Dobson. Patterns of macroparasite ag- gregation in wildlife host p opulations. Par asitolo gy , 117(6):597–608, 1998. doi:10.1017/S0031182098003448. [48] G. Stanek. P andora’s b ox: P athogens in Ixo des ricinus ticks in central Europ e [B ¨ uc hse der Pandora: Krankheitserreger in Ixo des ricinus -Zeck en in Mitteleuropa]. Wiener Klinische Wo chenschrift , 121(21-22):673–683, 2009. doi:10.1007/s00508- 009-1281-9. [49] A. Stradiotto, F. Cagnacci, R. Delaha y , S. Tioli, L. Nieder, and A. Rizzoli. Spatial organization of the yello w-nec ked mouse: Effects of density and resource av ailability . Journal of Mammalo gy , 90(3):704–714, 2009. doi:10.1644/08-MAMM-A-120R1.1. [50] D. Sumilo, L. Asokliene, A. Bormane, V. V asilenko, I. Golovljo v a, and S.E. Ran- dolph. Climate change cannot explain the upsurge of tic k-b orne encephalitis in the baltics. PL oS ONE , 2(6):e500, 2007. doi:10.1371/journal.p one.0000500. [51] Q.H. V uong. Likelihoo d ratio tests for mo del selection and non-nested h yp otheses. Ec onometric a , 57(2):307–333, 1989. 18 [52] Michael S. Williams and Eric D. Eb el. Metho ds for fitting the p oisson-lognormal distribution to microbial testing data. F o o d Contr ol , 27(1):73–80, 2012. ISSN 0956- 7135. doi:http://dx.doi.org/10.1016/j.foo dcon t.2012.03.007. [53] K. Wilson, Bj¨ ornstad. O.N., A.P . Dobson, S. Merler, G. Pogl ay en, S.E. Randolph, A.F. Read, and A Sk orping. Heterogeneities in macroparasite infections: patterns and pro cesses. In P . J. Hudson, A. Rizzoli, B.T. Grenfell, H. Heesterb eek, and A. P . Dobson, editors, The Ec olo gy of Wild life Dise ases , pages 6–44. Oxford Universit y Press, 2002. [54] W.L. Wilson, W.I. Mon tgomery , and R.W. Elwoo d. P opulation regulation in the w o o d mouse Ap o demus sylvaticus (L.). Mammal R eview , 23(2):73–92, 1993. [55] M.E.J. W o olhouse, C. Dye, J.F. Etard, T. Smith, J.D. Charlw o o d, G.P . Garnett, P . Hagan, J.L.K. Hii, P .D. Ndhlovu, R.J. Quinnell, C.H. W atts, S.K. Chandi- w ana, and R.M. Anderson. Heterogeneities in the transmission of infectious agen ts: Implications for the design of control programs. Pr o c e e dings of the Na- tional A c ademy of Scienc es of the Unite d States of Americ a , 94(1):338–342, 1997. doi:10.1073/pnas.94.1.338. [56] Z.-Y. Yin, H. Ren, Q.-M. Zhang, S.-L. Peng, Q.-F. Guo, and G.-Y. Zhou. Sp ecies abundance in a forest communit y in south chi na: A case of p oisson log- normal distribution. Journal of Inte gr ative Plant Biolo gy , 47(7):801–810, 2005. doi:10.1111/j.1744-7909.2005.00095.x. [57] Jerrold H. Zar. Biostatistic al Analysis . Prentice Hall, 2010. 19 Figures 2000 2001 2002 2003 2004 2005 2006 2007 2008 0 20 40 A.flavic ol lis A C E F B D G H I Figure 1: T emp oral v ariation of A.flavic ol lis mice abundance recorded in dif- feren t grids (lab elled in differen t colours from A to I). 20 1 θ J F M A M J J A S O N D Figure 2: Detection of seasonal abundance time-windo ws. The time series of captured mice has b een in terp olated b y a quadratic p olynomial curve. By normalising the obtained parab ola to unity and setting a threshold θ (= 0 . 5 in the example), we identify mice captured in high abundance season, those ab o v e the threshold θ (triangles), and mice captured in low abundance p erio d, those b elo w the threshold (circles). 21 0 0 . 1 0 . 2 0 . 3 0 . 4 0 20 40 60 NB PL P oiLN k PL min k NB min k PoiLN min k min KS 0 0 . 1 0 . 3 0 . 6 0 20 40 60 NB PL P oiLN k min GOF p -v alue 0 10 − 10 − 20 − 30 0 20 40 60 PL vs P oiLN PL vs NB k min LLR+V uong Figure 3: Comparison among fittings of distributions of ticks p er host with differen t functions. Left: Kolmogoro v-Smirnov statistic b et ween subsets of data ab o v e k min and the fitting mo dels on these subsets. V ertical dotted lines represent the optimum v alue of k min for different mo dels (NB: magenta; P oiLN: green; PL: cy an). F or the NB and P oiLN mo dels the optimum is observ ed for k NB min = k PoiLN min = 1, i.e. on the entire data set, while for the PL mo del the optimum is reac hed for k PL min = 38. Cen ter: goo dness-of-fit p - v alue of fitting mo dels on data larger than or equal to k min . As suggested b y Clauset and collab orators [13] for p -v alue greater than 0 . 1 (horizontal line) the fitting mo del is a go od description of the data. F or NB the GOF is lo w ( p < 10 − 3 ), suggesting the inappropriateness of the NB mo del in describing the data. The GOF of the PoiLN indicates that the mo del is appropriate only for large v alue of k min , th us simultaneously with large v alues of KS and therefore p oin ting out the lo w p erformance of the mo del. The PL fits should not b e rejected for v alues of k min larger than 35 concurrently with the low est v alue of KS. Right: Log-likelihoo d Ratio (LLR) test with V uong’s sign in terpretation. Negativ e (p ositiv e) v alues suggest the alternative mo del NB (red) or P oiLN (blue) distributions are (are not) fa voured in describing v alues larger than k min when compared to PL. The horizontal line shows the sign threshold. F ull marks sho w statistically significant tests ( p < 0 . 05) while empt y marks refer to non significan t tests ( p > 0 . 05). 22 1 10 100 1 10 − 1 10 − 2 10 − 3 10 − 4 P ( k ≥ x ) x (tic ks p er mouse) PL fit NB fit real-data P oiLN fit Figure 4: Complemen tary cum ulative functions of n umber of ticks p er host (real-data) with the b est p o wer-la w (PL), negativ e binomial (NB), and Poisson LogNormal (PoiLN) fit. The PL fitting mo del shows high pro ximity to the tail of the real data distribution while the NB and the PoiLN fits appropriately describ e the initial part of the distribution they describ e the tail improp erly . p 0 . 05 0 . 1 IN OUT PL PL α 2 3 4 IN OUT NB NB r 0 . 9 1 . 1 1 . 3 IN OUT k min 20 30 40 IN OUT Figure 5: Estimated parameters of differen t distributions (NB on left and PL on right) obtained inside (blue) and outside (red) of the mice p eak abundance time windo w. Time windows are defined by θ = 0 . 4 , 0 . 5 , 0 . 6 (from left to right for each subsets). V ertical bars indicate b est mo del fits (cen tral horizontal lines) with their uncertain ties that are 95% confidence in- terv al for NB mo dels while standard deviations for PL mo dels. 23 0 . 0 0 . 2 0 . 4 0 . 6 0 . 8 1 . 0 β 0 . 0 0 . 2 0 . 4 0 . 6 0 . 8 1 . 0 π ( t max ) f = 10% f = 5% f = 2% Figure 6: Median (line), interquartile (dark er area) and 95% confidence in- terv als (ligh ter area) of the final prev alence as a function of the transmission probability , for different v alues of f (= 2% , 5% , 10%), fraction of nymphs among ticks on a mouse, and by describing the ticks aggregation with the empirical distribution. Other parameters are N h = 10 4 , t max = 10 3 . 24 0 . 0 0 . 2 0 . 4 0 . 6 0 . 8 1 . 0 β 0 . 0 0 . 2 0 . 4 0 . 6 0 . 8 1 . 0 π ( t max ) P oiLN PL real NB Figure 7: Median (line), in terquartile (dark er area) and 95% confidence inter- v als (lighter area) of the final prev alence as a function of the trans- mission probability , for differen t fitting distributions (PL, NB and PoiLN). Other parameters are N h = 10 4 , t max = 10 3 , and f = 2%. 25 0 . 0 0 . 2 0 . 4 0 . 6 0 . 8 1 . 0 β 0 . 0 0 . 2 0 . 4 0 . 6 0 . 8 1 . 0 π ( t max ) k min = 5 k min = 10 k min = 15 PL PL with temporal windo w Figure 8: Median of the final prev alence as a function of the transmission probabilit y . A PL distribution of vectors-per-host has b een considered in all scenarios. Simulations that consider different aggregation b ehaviours accord- ing to the temporal window of mice abundance (red) are compared with others with a fixed distribution (blue). Other parameters are N h = 10 4 , t max = 10 3 , and f = 10%. T ables T able 1: Basic descriptiv e statistics for empirical data. 2000 2001 2002 2003 2004 2005 2006 2007 2008 # of grids 9 8 8 4 4 4 4 4 4 # of sessions 16 15 15 13 10 15 15 16 14 # of mice 1207 356 434 137 187 854 327 897 323 sum of feeding ticks 14376 7073 6550 2426 3063 7361 4077 5821 4685 median of ticks p er ro den t 9 14 11 14 11 6 8 4 10 ranges of ticks p er ro den t (1 , 103) (1 , 102) (1 , 78) (3 , 88) (1 , 95) (1 , 111) (1 , 93) (1 , 85) (1 , 77) nymphs fraction 0 . 5% 2 . 2% 2 . 7% 1 . 1% 3 . 0% 0 . 8% 2 . 0% 0 . 7% 1 . 0% Num b er of trapping grids, trapping sessions, total n umber of A. flavic ol lis captures for differen t years, sum of feeding ticks, median and ranges of the num b er of ticks p er ro den t, and mean n umber of nymphs fraction among feeding ticks. 26 T able 2: Nymphs to total tic ks ratio for observed feeding tic ks on mice. 2000 2001 2002 2003 2004 2005 2006 2007 2008 A 1 . 1 1 . 0 4 . 0 0 . 3 2 . 0 1 . 1 1 . 8 0 . 6 1 . 6 B 0 . 5 1 . 6 0 . 7 1 . 0 4 . 6 0 . 6 1 . 8 0 . 7 0 . 1 C 0 . 4 1 . 2 2 . 8 1 . 3 1 . 7 0 . 4 3 . 1 1 . 2 1 . 4 D 0 . 2 1 . 8 3 . 6 1 . 0 3 . 9 0 . 9 1 . 4 0 . 5 0 . 5 E 0 . 4 4 . 4 1 . 6 - - - - - - F 0 . 9 2 . 8 2 . 2 - - - - - - G 0 . 3 - - - - - - - - H 1 . 1 1 . 5 3 . 0 - - - - - - I 0 . 7 10 . 4 3 . 3 - - - - - - P ercentage of feeding nymphs on the total feeding ticks observed on mice in different y ears (2000-2008) and grids ( A - I ). SM.1 Sensitivit y Analysis of Distribution P a rameters on Epidemic Spreading In this section we explored the effect of parameters of tick burden distributions not at the b est fit on the epidemic spreading. Results suggested that the larger the heterogeneity caused b y parameters, the low er the epidemic threshold. 27 0 . 0 0 . 2 0 . 4 0 . 6 0 . 8 1 . 0 β 0 . 0 0 . 1 0 . 2 0 . 3 0 . 4 0 . 5 0 . 6 π ( t max ) PL with α − 1 PL with α PL with α + 1 real Figure 9: Median (lines), interquartile (dark er areas) and 95% confidence in- terv als (lighter areas) of the final prev alence as a function of the transmission probability . Ticks burdens are describ ed by PL distribution with different exp onen ts (see legend). In particular, we explore the sensitivit y of this function to v ariations in α , which represen ts the b est fit parameter on the empirical data. The long-term prev alence obtained with tick burdens sampled from the empirical distribution is also plotted as b enc hmark. 28 0 . 0 0 . 2 0 . 4 0 . 6 0 . 8 1 . 0 β 0 . 0 0 . 1 0 . 2 0 . 3 0 . 4 0 . 5 0 . 6 π ( t max ) NB with r − 0 . 5 and p − 0 . 05 NB with r and p − 0 . 05 NB with r + 0 . 5 and p − 0 . 05 NB with r − 0 . 5 and p NB with r and p NB with r + 0 . 5 and p NB with r − 0 . 5 and p + 0 . 05 NB with r and p + 0 . 05 NB with r + 0 . 5 and p + 0 . 05 real Figure 10: Median (line), in terquartile (dark er areas) and 95% confidence in- terv als (ligh ter areas) of the final prev alence as a function of the transmission probability . Tic ks burdens are describ ed b y NB distribution with differen t parameters (see legend). In particular, we explore the sensi- tivit y of this function to v ariations in r and p , which represen t the b est fit parameters on the empirical data. The final prev alence obtained with tick burdens sampled from the empirical distribution is also plotted as b enc hmark. 29 0 . 0 0 . 2 0 . 4 0 . 6 0 . 8 1 . 0 β 0 . 0 0 . 1 0 . 2 0 . 3 0 . 4 0 . 5 0 . 6 π ( t max ) P oiLN with µ − 0 . 5 and σ − 0 . 2 P oiLN with µ and σ − 0 . 2 P oiLN with µ + 0 . 5 and σ − 0 . 2 P oiLN with µ − 0 . 5 and σ P oiLN with µ and σ P oiLN with µ + 0 . 5 and σ P oiLN with µ − 0 . 5 and σ + 0 . 2 P oiLN with µ and σ + 0 . 2 P oiLN with µ + 0 . 5 and σ + 0 . 2 real Figure 11: Median (line), in terquartile (dark er areas) and 95% confidence in- terv als (ligh ter areas) of the final prev alence as a function of the transmission probability . Ticks burdens are describ ed b y P oiLN distri- bution with differen t couple of parameters. In particular, we explore the sensitivit y of this curve to v ariations of µ and σ , the b est fit parameters on the empirical data. The final prev alence obtained with tick burdens sampled from the empirical distribution is also plotted as b enc hmark. 30 SM.2 Sensitivit y Analysis of f , F raction of Nymphs among Ticks In this section we explored the effect of differen t fractions, f , of n ymphs ov erall the total n umber of tic ks on epidemic spreading. By exploring f = 2% in main text, f = 5% in Figure 12 and f = 10% in Figure 13 w e conclude that the larger the f , and the larger the probability of inv asion of the pathogen. Moreov er, it is worth to stress out that for differen t v alues of f epidemic curves qualitativ ely do not change and in particular the order of epidemic thresholds is main tained. 0 . 0 0 . 2 0 . 4 0 . 6 0 . 8 1 . 0 β 0 . 0 0 . 2 0 . 4 0 . 6 0 . 8 1 . 0 π ( t max ) P oiLN PL real NB Figure 12: Median (lines), interquartile (darker areas) and 95% confidence in- terv als (ligh ter areas) of the final prev alence as a function of the transmission probabilit y . f , fraction of n ymphs among ticks is fixed to 5%. 31 0 . 0 0 . 2 0 . 4 0 . 6 0 . 8 1 . 0 β 0 . 0 0 . 2 0 . 4 0 . 6 0 . 8 1 . 0 π ( t max ) P oiLN PL real NB Figure 13: Median (lines), interquartile (darker areas) and 95% confidence in- terv als (ligh ter areas) of the final prev alence as a function of the transmission probabilit y . f , fraction of n ymphs among ticks is fixed to 10%. 32