Open Access

Impact of biodiversity and seasonality on Lyme-pathogen transmission

Theoretical Biology and Medical Modelling201411:50

https://doi.org/10.1186/1742-4682-11-50

Received: 1 August 2014

Accepted: 22 October 2014

Published: 28 November 2014

Abstract

Lyme disease imposes increasing global public health challenges. To better understand the joint effects of seasonal temperature variation and host community composition on the pathogen transmission, a stage-structured periodic model is proposed by integrating seasonal tick development and activity, multiple host species and complex pathogen transmission routes between ticks and reservoirs. Two thresholds, one for tick population dynamics and the other for Lyme-pathogen transmission dynamics, are identified and shown to fully classify the long-term outcomes of the tick invasion and disease persistence. Seeding with the realistic parameters, the tick reproduction threshold and Lyme disease spread threshold are estimated to illustrate the joint effects of the climate change and host community diversity on the pattern of Lyme disease risk. It is shown that climate warming can amplify the disease risk and slightly change the seasonality of disease risk. Both the “dilution effect” and “amplification effect” are observed by feeding the model with different possible alternative hosts. Therefore, the relationship between the host community biodiversity and disease risk varies, calling for more accurate measurements on the local environment, both biotic and abiotic such as the temperature and the host community composition.

Keywords

Seasonal tick populationLyme diseaseHost diversityClimateThresholdDilution effectAmplification effect

Introduction

Lyme disease is acknowledged as a common infectious disease for the most of the world, especially in Europe and North America. The disease is caused by a bacterium called Borrelia burgdorferi, transmitted by ticks, especially Ixodes scapularis[1, 2]. It affects both humans and animals, with more than 30,000 cases reported annually in the United States alone [3]. The pathogen transmission involves three ecological and epidemiological processes: nymphal ticks infected in the previous year appear first; these ticks then transmit the pathogen to their susceptible vertebrate hosts during a feeding period; the next generation larvae acquire infection by sucking recently infected hosts’ blood and these larvae develop into nymphs in the next year to complete the transmission cycle.

Understanding the factors that regulate the abundance and distribution of the Lyme-pathogen is crucial for the effective control and prevention of the disease. Host diversity and temperature variation have direct influence on Lyme transmission patterns [4]. The tick vectors need to complete the transition of four life stages of metamorphosis (eggs, larvae, nymphs and adults) and each postegg stage requires a blood meal from a wide range of host species, and every host species has a specialized reservoir competence, namely ability to carry and transmit the pathogen [5]. Moreover, weather conditions (temperature, rainfall, humidity, for example) are known to affect the reproduction, development, behavior, and population dynamics of the arthropod vectors [69], thereby the spread of the Lyme-pathogen in vectors. In particular, the temperature is regarded as an important factor affecting the tick development and tick biting activity, which gives rise to tick seasonal dynamics [1, 10, 11]. In summary, host diversity [1216], stage structure of ticks [1, 13, 1722] and climate effects [1, 10, 11, 13, 20, 23] are considered to be crucial for the persistence of Lyme infection. Therefore modeling Lyme-pathogen transmission with multiple tick life stages, tick seasonality and host community composition is pivotal in understanding the pathogen transmission.

There have been a range of tick-borne disease modeling efforts dedicating to different aspects of Lyme disease transmission: the basic Lyme transmission ecology [24, 25], effect of different hosts and their densities on the persistence of tick-borne diseases [15, 16, 26], threshold dynamics for disease infection [27, 28], seasonal tick population dynamics and disease transmission [2, 10, 29], climatic effects [1, 2, 22], spatial heterogeneity [18, 30, 31], among others. These previous studies promote our understanding on the transmission mechanisms and designing effective prevention and control measures. In this paper, we develop a modeling framework incorporating the impact of multiple tick life stages, tick seasonality and host diversity on the Lyme disease transmission cycle. We follow the generic model proposed by Randolph and Rogers [21], and divide the vector population into 7 stages with 12 subclasses, as illustrated in Figure 1. This generic model is able to account for the following key features: (i) temperature-dependent/temperature-independent development rates; (ii) temperature-dependent host seeking rates; (iii) density-dependent mortalities, caused by the hosts’ responses during the feeding period; (iv) constant mortalities of off-host development stages. The proposed model below is different from these existing models by incorporating all aforementioned aspects of Lyme disease transmission in a single framework, and as such this framework permits us to analytically define the thresholds of tick population dynamics and pathogen transmission dynamics under seasonal temperature variation, and establish the relationship of these thresholds to the tick establishment and pathogen persistence.
Figure 1

Schematic diagram for the Lyme disease transmission. To describe the tick development and biting activities, the tick population is divided into 7 stages, stratified further as the uninfected or infected epidemiological classes for postegg stages. Immature ticks can feed on two host species, the mice (H1) and an alternative host (H2), while adult ticks are assumed to feed only on deer in this study.

The remaining parts of the paper are organized as follows. In the spirit of striking a delicate balance between the feasibility for the recognized mathematical analysis and the necessity for capturing the key ecological/epidemiological reality, a stage-structured deterministic model is formulated in the next section. Moreover, two thresholds, one for the tick population dynamics and the other for the Lyme-pathogen transmission dynamics are derived and shown to be pivotal in determining the tick population establishment and disease invasion. The model is parameterized in section ‘Model parametrization’. Then the question concerning whether the climate change and an additional host species can amplify/dilute disease prevalence and change the seasonality of disease risk will be addressed through model simulation in section ‘Results’. A discussion section concludes the paper.

Mathematical model

Model formulation

In line with the complex physiological process of Ixodes scapularis, we divide them into four stages: eggs (E), larvae (L), nymphs (N) and adults (A). Each postegg stage is further divided into two groups: questing (Q) and feeding (F) according to their behavior on or off hosts. Moreover, in terms of their infection status, each group is stratified into two subgroups: susceptible (S) and infected (I). All variable notations are self-explained as summarized in Table 1. For instance, L FS represents the subgroup of susceptible feeding larvae.
Table 1

Variable explanations used in the model (1)

Variable

Meaning

E

the number of eggs

L Q

the number of questing larvae

L FS

the number of susceptible feeding larvae

L FI

the number of infected feeding larvae

N QS

the number of susceptible questing nymphs

N QI

the number of infected questing nymphs

N FS

the number of susceptible feeding nymphs

N FI

the number of infected feeding nymphs

A QS

the number of susceptible questing adults

A QI

the number of infected questing adults

A FS

the number of susceptible feeding adults

A FI

the number of infected feeding adults

H 1I

the number of infected white-footed mice

H 2I

the number of infected alternative hosts

We assume that the host community of the tick population contains three species groups: (i) the white-footed mice H1 (mainly Peromyscus leucopus) with the mortality rate μ H 1 , which is widely known as a primary food provider of immature I. scapularis ticks and a key reservoir competent host of B. burgdorferi reflecting the strong ability to be infected with the pathogen and to transmit the pathogen to its vector; (ii) the white-tailed deer D (mainly Odocoileus virginianus), which is believed to be the paramount food provider for adults and in-transmissible for the spread of Lyme-pathogen [32]; and (iii) an alternative host H2 with mortality rate μ H 2 such as the eastern chipmunk, the Virginia opossum and the western fence lizard, which is used to study the impact of host community composition on the Lyme disease risk. For the sake of simplicity, we further assume that the total number of each host species (susceptible plus infected) in an isolated habitat is constant. However the number of infected hosts can vary with time, denote by H1I and H2I, respectively. The Lyme-pathogen transmission cycle between the hosts and multi-stage tick population is presented in the diagram of Figure 1.

In the host-pathogen-tick transmission cycle, larvae and nymphs will bite their host species, however their biting preference to different host species may be different. In order to identify the difference, we use the coefficients, p1(p2), to describe larval (nymphal) ticks biting bias on their hosts [33, 34]. Specifically, p1>1(p2>1) indicates one host H2 can attract more larval (nymphal) bites than one host H1 and vice versa when 0<p1<1(0<p2<1). Using the method described in [35], F L ( t ) H 1 H 1 + p 1 H 2 H 1 I ( t ) H 1 is the average rate at which a susceptible questing larva finds and attaches successfully onto the infected mice, where F L (t) is the feeding rate of larvae, and then β H 1 L F L ( t ) H 1 H 1 + p 1 H 2 H 1 I ( t ) H 1 is the average infection rate at which a susceptible larva gets infected from mice, where β H 1 L is the pathogen transmission probability per bite from infectious mice H1 to susceptible larvae. Using the same idea, the infection rate of larvae from the infected alternative host H2 can be accounted. Therefore, the larval infection rate is given by
β H 1 L F L ( t ) H 1 H 1 + p 1 H 2 H 1 I ( t ) H 1 L Q ( t ) + β H 2 L F L ( t ) p 1 H 2 H 1 + p 1 H 2 H 2 I ( t ) H 2 L Q ( t ) = β H 1 L H 1 I ( t ) H 1 + p 1 H 2 + β H 2 L p 1 H 2 I ( t ) H 1 + p 1 H 2 F L ( t ) L Q ( t ) .
Similarly, the nymphal infection rate which comes from the contact of questing susceptible nymphs and infectious hosts is given by
β H 1 N F N ( t ) H 1 H 1 + p 2 H 2 H 1 I ( t ) H 1 N QS ( t ) + β H 2 N F N ( t ) p 2 H 2 H 1 + p 2 H 2 H 2 I ( t ) H 2 N QS ( t ) = β H 1 N H 1 I ( t ) H 1 + p 2 H 2 + β H 2 N p 2 H 2 I ( t ) H 1 + p 2 H 2 F N ( t ) N QS ( t ) .
The susceptible hosts can get infected when they are bitten by infected questing nymphs. The conservation of bites requires that the numbers of bites made by ticks and received by hosts should be the same. The disease incidence rate for mice is therefore given by
F N ( t ) β N H 1 ( N QI ( t ) + N QS ( t ) ) N QI ( t ) N QI ( t ) + N QS ( t ) H 1 H 1 + p 2 H 2 H 1 H 1 I ( t ) H 1 = F N ( t ) β N H 1 N QI ( t ) H 1 H 1 I ( t ) H 1 + p 2 H 2 .
Similarly, the alternative host is infected by the infectious nymphal biting at a rate
F N ( t ) β N H 2 N QI ( t ) p 2 ( H 2 H 2 I ( t ) ) H 1 + p 2 H 2 .
Therefore, the disease transmission process between ticks and their hosts can be described by the following system:
dE dt = b ( t ) ( A FS ( t ) + A FI ( t ) ) μ E E ( t ) d E ( t ) E ( t ) , d L Q dt = d E ( t ) E ( t ) μ QL L Q ( t ) F L ( t ) L Q ( t ) , d L FS dt = 1 β H 1 L H 1 I ( t ) H 1 + p 1 H 2 + β H 2 L p 1 H 2 I ( t ) H 1 + p 1 H 2 F L ( t ) L Q ( t ) μ FL L FS ( t ) D L ( L FS ( t ) + L FI ( t ) ) L FS ( t ) d L ( t ) L FS ( t ) , d L FI dt = β H 1 L H 1 I ( t ) H 1 + p 1 H 2 + β H 2 L p 1 H 2 I ( t ) H 1 + p 1 H 2 F L ( t ) L Q ( t ) μ FL L FI ( t ) D L ( L FS ( t ) + L FI ( t ) ) L FI ( t ) d L ( t ) L FI ( t ) , d N QS dt = d L ( t ) L FS ( t ) μ QN N QS ( t ) F N ( t ) N QS ( t ) , d N QI dt = d L ( t ) L FI ( t ) μ QN N QI ( t ) F N ( t ) N QI ( t ) , d N FS dt = 1 β H 1 N H 1 I ( t ) H 1 + p 2 H 2 + β H 2 N p 2 H 2 I ( t ) H 1 + p 2 H 2 F N ( t ) N QS ( t ) μ FN N FS ( t ) D N ( N FS ( t ) + N FI ( t ) ) N FS ( t ) d N ( t ) N FS ( t ) , d N FI dt = F N ( t ) N QI ( t ) + β H 1 N H 1 I ( t ) H 1 + p 2 H 2 + β H 2 N p 2 H 2 I ( t ) H 1 + p 2 H 2 F N ( t ) N QS ( t ) μ FN N FI ( t ) D N ( N FS ( t ) + N FI ( t ) ) N FI ( t ) d N ( t ) N FI ( t ) , d A QS dt = d N ( t ) N FS ( t ) μ QA A QS ( t ) F A ( t ) A QS ( t ) , d A QI dt = d N ( t ) N FI ( t ) μ QA A QI ( t ) F A ( t ) A QI ( t ) , d A FS dt = F A ( t ) A QS ( t ) μ FA A FS ( t ) D A ( A FS ( t ) + A FI ( t ) ) A FS ( t ) , d A FI dt = F A ( t ) A QI ( t ) μ FA A FI ( t ) D A ( A FS ( t ) + A FI ( t ) ) A FI ( t ) , d H 1 I dt = F N ( t ) β N H 1 N QI ( t ) H 1 H 1 I ( t ) H 1 + p 2 H 2 μ H 1 H 1 I ( t ) , d H 2 I dt = F N ( t ) β N H 2 N QI ( t ) p 2 ( H 2 H 2 I ( t ) ) H 1 + p 2 H 2 μ H 2 H 2 I ( t ) .
(1)
We assume all the coefficients in the system are nonnegative and the time-dependent coefficients are τ-periodic with period τ=365 days. The detailed parameter definitions and sample values of these parameters are represented in Table 2. These time-dependent parameters will be estimated in subsection ‘Time-dependent parameters’ below.
Table 2

Definitions and corresponding values of the model parameter with the daily timescale

Parameter

Meaning

(Value, [reference]) or estimation

μ E

mortality rate of eggs

(0.0025, [1])

μ QL

mortality rate of questing larvae

(0.006, [1])

μ QN

mortality rate of questing nymphs

(0.006, [1])

μ QA

mortality rate of questing adults

(0.006, [1])

μ FL

natural mortality rate of feeding larvae

(0.038, A)

μ FN

natural mortality rate of feeding nymphs

(0.028, A)

μ FA

natural mortality rate of feeding adults

(0.018, A)

H 1

the number of white-footed mice

(200, [1])

β H 1 L

transmission probability from H1 to larvae

(0.6, [13])

β N H 1

transmission probability from nymphs to H1

(1, [13])

μ H 1

death rate of the white-footed mice

(0.012, [13])

H 2

the number of alternative host H2

(variable)

β H 2 L

transmission probability from H2 to larvae

(variable, [36])

β N H 2

transmission probability from nymphs to H2

(variable, [36])

μ H 2

death rate of the alternative host H2

(variable)

D

the number of deer

(20, [1])

p

the maximum number of eggs produced

(3000, [1])

p 1

larval biting bias for host H2

(variable, [37])

p 2

nymphal biting bias for host H2

(variable, [37])

b(t)

birth rate of eggs produced

(see subsection ‘Time-dependent parameters’)

d E (t)

development rate of eggs

(see subsection ‘Time-dependent parameters’)

d L (t)

development rate of larvae

(see subsection ‘Time-dependent parameters’)

d N (t)

development rate of nymphs

(see subsection ‘Time-dependent parameters’)

F L (t)

feeding rate of larvae

(see subsection ‘Time-dependent parameters’)

F N (t)

feeding rate of nymphs

(see subsection ‘Time-dependent parameters’)

F A (t)

feeding rate of adults

(see subsection ‘Time-dependent parameters’)

D L

density-dependent mortality rate of feeding larvae

( 0.001084 H 1 + p 1 H 2 , E)

D N

density-dependent mortality rate of feeding nymphs

( 0.001084 H 1 + p 2 H 2 , E)

D A

density-dependent mortality rate of feeding adults

( 0.001084 D , E)

Where E: estimation based on [38] and A: assumption.

Dynamics analysis

Positivity and boundedness of solutions

Our first task is to show that the mathematical model (1) is biologically meaningful. To do this, we first establish the following theorem to ensure that all solutions through nonnegative initial values remain nonnegative and bounded. You may refer Appendix 1 for the proof.

Theorem 2.1.
For each initial value x 0 X : = { x + 14 : H 1 I H 1 , H 2 I H 2 } , system (1) has a unique and bounded solution x(t,x0). Moreover, the solution x(t,x0)remains in X for any t≥0. Here, x + 14 denotes a generic point with components
( E , L Q , L FS , L FI , N QS , N QI , N FS , N FI , A QS , A QI , A FS , A FI , H 1 I , H 2 I ) .
Using change of variables L F =L FS +L FI , N Q =N QS +N QI , N F =N FS +N FI , A Q =A QS +A QI and A F =A FS +A FI , system (1) reduces to
dE dt = b ( t ) A F ( t ) ( μ E + d E ( t ) ) E ( t ) , d L Q dt = d E ( t ) E ( t ) ( μ QL + F L ( t ) ) L Q ( t ) , d L F dt = F L ( t ) L Q ( t ) D L L F 2 ( t ) ( μ FL + d L ( t ) ) L F ( t ) , d N Q dt = d L ( t ) L F ( t ) ( μ QN + F N ( t ) ) N Q ( t ) , d N F dt = F N ( t ) N Q ( t ) D N N F 2 ( t ) ( μ FN + d N ( t ) ) N F ( t ) , d A Q dt = d N ( t ) N F ( t ) ( μ QA + F A ( t ) ) A Q ( t ) , d A F dt = F A ( t ) A Q ( t ) μ FA A F ( t ) D A A F 2 ( t ) , d L FI dt = β H 1 L H 1 I ( t ) H 1 + p 1 H 2 + β H 2 L p 1 H 2 I ( t ) H 1 + p 1 H 2 F L ( t ) L Q ( t ) D L L F ( t ) L FI ( t ) ( μ FL + d L ( t ) ) L FI ( t ) , d N QI dt = d L ( t ) L FI ( t ) ( μ QN + F N ( t ) ) N QI ( t ) , d H 1 I dt = F N ( t ) β N H 1 N QI ( t ) H 1 H 1 I ( t ) H 1 + p 2 H 2 μ H 1 H 1 I ( t ) , d H 2 I dt = F N ( t ) β N H 2 N QI ( t ) p 2 ( H 2 H 2 I ( t ) ) H 1 + p 2 H 2 μ H 2 H 2 I ( t ) .
(2)

Note that we have other three equations for infected feeding nymphs (N FI ), questing adults (A QI ) and feeding adults (A FI ), which is decoupled from the above system. Biologically, we pay attention to the population size of infected questing nymphs whose bites are the main courses of human Lyme disease. We thereby focus on system (2) in the remaining of the paper.

The tick population dynamics

We firstly consider the following stage-structured system for the tick population growth decoupled from system (2):
dE dt = b ( t ) A F ( t ) ( μ E + d E ( t ) ) E ( t ) , d L Q dt = d E ( t ) E ( t ) ( μ QL + F L ( t ) ) L Q ( t ) , d L F dt = F L ( t ) L Q ( t ) D L L F 2 ( t ) ( μ FL + d L ( t ) ) L F ( t ) , d N Q dt = d L ( t ) L F ( t ) ( μ QN + F N ( t ) ) N Q ( t ) , d N F dt = F N ( t ) N Q ( t ) D N N F 2 ( t ) ( μ FN + d N ( t ) ) N F ( t ) , d A Q dt = d N ( t ) N F ( t ) ( μ QA + F A ( t ) ) A Q ( t ) , d A F dt = F A ( t ) A Q ( t ) μ FA A F ( t ) D A A F 2 ( t ) .
(3)
Linearization of system (3) at zero leads to the following linear system
dE dt = b ( t ) A F ( t ) ( μ E + d E ( t ) ) E ( t ) , d L Q dt = d E ( t ) E ( t ) ( μ QL + F L ( t ) ) L Q ( t ) , d L F dt = F L ( t ) L Q ( t ) ( μ FL + d L ( t ) ) L F ( t ) , d N Q dt = d L ( t ) L F ( t ) ( μ QN + F N ( t ) ) N Q ( t ) , d N F dt = F N ( t ) N Q ( t ) ( μ FN + d N ( t ) ) N F ( t ) , d A Q dt = d N ( t ) N F ( t ) ( μ QA + F A ( t ) ) A Q ( t ) , d A F dt = F A ( t ) A Q ( t ) μ FA A F ( t ) .
(4)
Let F(t)=(f ij (t))7×7, where f1,7(t)=b(t) and fi,j(t)=0 if (i,j)≠(1,7), and V(t) =
μ E + d E ( t ) 0 0 0 0 0 0 d E ( t ) μ QL + F L ( t ) 0 0 0 0 0 0 F L ( t ) μ FL + d L ( t ) 0 0 0 0 0 0 d L ( t ) μ QN + F N ( t ) 0 0 0 0 0 0 F N ( t ) μ FN + d N ( t ) 0 0 0 0 0 0 d N ( t ) μ QA + F A ( t ) 0 0 0 0 0 0 F A ( t ) μ FA .
Then we can rewrite (4) as
dx ( t ) dt = ( F ( t ) V ( t ) ) x ( t ) ,
where a vector x(t)=(E(t),L Q (t),L F (t),N Q (t),N F (t),A Q (t),A F (t)) T . Assume Y(t,s), ts, is the evolution operator of the linear periodic system dy dt = V ( t ) y . That is, for each s , the 7×7 matrix Y(t,s) satisfies
d dt Y ( t , s ) = V ( t ) Y ( t , s ) t s , Y ( s , s ) = I ,
where I is the 7×7 identity matrix. Let C τ be the Banach space of all τ-periodic functions from to 7 , equipped with the maximum norm. Suppose ϕC τ is the initial distribution of tick individuals in this periodic environment. Then F(s)ϕ(s) is the rate of new ticks produced by the initial ticks who were introduced at time s, and Y(t,s)F(s)ϕ(s) represents the distribution of those ticks who were newly produced at time s and remain alive at time t for ts. Hence,
ψ ( t ) = t Y ( t , s ) F ( s ) ϕ ( s ) ds = 0 Y ( t , t a ) F ( t a ) ϕ ( t a ) da

is the distribution of accumulative ticks at time t produced by all those ticks ϕ(s) introduced at the previous time.

Following ideas proposed in [39, 40], we define a next generation operator G:C τ C τ by
( G ϕ ) ( t ) = 0 Y ( t , t a ) F ( t a ) ϕ ( t a ) da t , ϕ C τ .

Then the spectral radius of G is defined as R v : = ρ ( G ) . In what follows, we call R v as a threshold for tick population dynamics.

Let Φ P (t) and ρ(Φ P (τ)) be the monodromy matrix of the linear τ-periodic system dz dt = P ( t ) z and the spectral radius of Φ P (τ), respectively. Then, from [40], Theorem 2.2, we conclude (i) R v = 1 if and only if ρ(ΦFV(τ))=1; (ii) R v > 1 if and only if ρ(ΦFV(τ))>1; (iii) R v < 1 if and only if ρ(ΦFV(τ))<1. We also know that the zero solution is locally asymptotically stable if R v < 1 , and unstable if R v > 1 .

Note that the Poincar e ́ map associated with system (3) is not strongly monotone since some coefficients are not strictly positive (remain zero in a nonempty interval). However, if we regard a τ-periodic system (3) as a 6τ-periodic system, we can show that the Poincar e ́ map with respect to the 6τ-periodic system is strongly monotone by using the same idea as in [41], Lemma 3.2. We then use [42], Theorem 2.3.4, to the Poincar e ́ map associated with system (3) to obtain the following result, with the proof in Appendix 2.

Theorem 2.2.
The following statements are valid:
  1. (i)

    If R v 1 , then zero is globally asymptotically stable for system (3) in + 7 ;

     
  2. (ii)
    If R v > 1 , then system (3) admits a unique τ-positive periodic solution
    ( E ( t ) , L Q ( t ) , L F ( t ) , N Q ( t ) , N F ( t ) , A Q ( t ) , A F ( t ) ) ,
     

and it is globally asymptotically stable for system (3) with initial values in + 7 { 0 } .

The global dynamics of the full model

If threshold for ticks R v > 1 , then there exists a positive periodic solution,
( E ( t ) , L Q ( t ) , L F ( t ) , N Q ( t ) , N F ( t ) , A Q ( t ) , A F ( t ) ) ,
for system (3) such that
lim t { ( E ( t ) , L Q ( t ) , L F ( t ) , N Q ( t ) , N F ( t ) , A Q ( t ) , A F ( t ) ) ( E ( t ) , L Q ( t ) , L F ( t ) , N Q ( t ) , N F ( t ) , A Q ( t ) , A F ( t ) ) } = 0 .
In this case, equations for the infected populations in system (2) give rise to the following limiting system:
d L FI dt = β H 1 L H 1 I ( t ) H 1 + p 1 H 2 + β H 2 L p 1 H 2 I ( t ) H 1 + p 1 H 2 F L ( t ) L Q ( t ) D L L F ( t ) L FI ( t ) ( d L ( t ) + μ FL ) L FI ( t ) , d N QI dt = d L ( t ) L FI ( t ) ( μ QN + F N ( t ) ) N QI ( t ) , d H 1 I dt = F N ( t ) β N H 1 N QI ( t ) H 1 H 1 I ( t ) H 1 + p 2 H 2 μ H 1 H 1 I ( t ) , d H 2 I dt = F N ( t ) β N H 2 N QI ( t ) p 2 ( H 2 H 2 I ( t ) ) H 1 + p 2 H 2 μ H 2 H 2 I ( t ) .
(5)
Following ideas of [39, 40], as proceed in the definition of R v in the previous section, we can define a threshold for the pathogen. To do this, we introduce
F ~ ( t ) = 0 0 β H 1 L F L ( t ) L Q ( t ) H 1 + p 1 H 2 p 1 β H 2 L F L ( t ) L Q ( t ) H 1 + p 1 H 2 0 0 0 0 0 F N ( t ) β N H 1 H 1 H 1 + p 2 H 2 0 0 0 p 2 F N ( t ) β N H 2 H 2 H 1 + p 2 H 2 0 0 ,
and
V ~ ( t ) = D L L F ( t ) + d L ( t ) + μ FL 0 0 0 d L ( t ) μ QN + F N ( t ) 0 0 0 0 μ H 1 0 0 0 0 μ H 2 .
Assume Y ~ (t,s), ts, is the evolution operator of the linear periodic system dy dt = V ~ ( t ) y . Let C ~ τ be the Banach space of all τ-periodic functions from to 4 , equipped with the maximum norm. Suppose ϕ C ~ τ is the initial distribution of infectious tick and host individuals in this periodic environment. Then F ~ ( s ) ϕ ( s ) is the rate of new infectious ticks and host individuals produced by the initial infectious ticks and hosts who were introduced at time s, and Y ~ ( t , s ) F ~ ( s ) ϕ ( s ) represents the distribution of those ticks who were newly produced at time s and remain alive at time t for ts. Hence,
ψ ~ ( t ) = t Y ~ ( t , s ) F ~ ( s ) ϕ ( s ) ds = 0 Y ~ ( t , t a ) F ~ ( t a ) ϕ ( t a ) da
is the distribution of accumulative infectious ticks and hosts at time t produced by all those infectious individuals ϕ(s) introduced at the previous time. Define the a next generation operator G ~ : C ~ τ C ~ τ by
( G ~ ( ϕ ) ) ( t ) = 0 Y ~ ( t , t a ) F ~ ( t a ) ϕ ( t a ) da t , ϕ C τ ~ .

It then follows from [39, 40] that the spectral radius of G ~ is define as R 0 : = ρ ( G ~ ) , and shows that it is a threshold of the Lyme-pathogen dynamics (5).

Using the same argument as in the proof of Theorem 2.2 (see also the proof of Lemma 2.3 in [43]), we have the following results:

Theorem 2.3.

(i) If R 0 1 , then zero is globally asymptotically stable for system (5) in + 4 ; (ii) If R 0 > 1 , then system (5) admits a unique positive periodic solution ( L FI ( t ) , N QI ( t ) , H 1 I ( t ) , H 2 I ( t ) ) and it is globally asymptotically stable for system (5).

Based on the aforementioned two thresholds, R v for ticks dynamics and R 0 for the pathogen dynamics, we can completely determine the global dynamics of the system (2). The detailed proof is shown in Appendix 3.

Theorem 2.4.

Let x(t,x0) be the solution of system (2) through x0. Then the following statements are valid:
  1. (i)

    If R v 1 , then zero is globally attractive for system (2);

     
  2. (ii)
    If R v > 1 and R 0 1 , then
    lim t { ( x 1 ( t ) , x 2 ( t ) , x 3 ( t ) , x 4 ( t ) , x 5 ( t ) , x 6 ( t ) , x 7 ( t ) ) ( E ( t ) , L Q ( t ) , L F ( t ) , N Q ( t ) , N F ( t ) , A Q ( t ) , A F ( t ) ) } = 0 ,

    and lim t x i ( t ) = 0 for i[8,11];

     
  3. (iii)

    If R v > 1 and R 0 > 1 , then there exists a positive periodic solution x (t), and this periodic solution is globally attractive for system (2) with respect to all positive solutions.

     

Summary of mathematical results

By incorporating the tick physiological development and multiple host species, we propose a seasonal deterministic stage-structured Lyme disease transmission model. The model turns out to be a periodic system of ordinary differential equations with high dimensions. As the pathogen has a negligible effect on population dynamics of the ticks and their hosts, the dynamics of the ticks is independent of the pathogen occurrence. This allows us to obtain an independent subsystem for the dynamics of the tick population. Taking the advantage of this observation and with the help of the developed theory for chain transitive sets, we are able to derive two results on global stability of the model system (2). Two biologically significant indices, the tick reproduction threshold R v and the Lyme disease invasion threshold R 0 are derived and shown to completely classify the long term outcomes of the tick and pathogen establishment.

Model parametrization

In this section, we present the estimation of the time-dependent parameters and other parameters related to the host species.

Alternative hosts species and their reservoir competence

To study the potential effect of of host community biodiversity on the risk of Lyme disease, three types of alternative host species are considered which are different from their reservoir competence, namely, the product of host infection probability bitten by infectious nymphs and larvae infection probability from infectious hosts [36]. The first type is considered as the one with high reservoir competence such as the short-tailed shrew, the marked shrew and the eastern chipmunk. The values of β H 2 L and β N H 2 are set as 0.569 and 0.971, respectively, as reported in [36]. The second type that we want to compare is the one with low reservoir competence, in which β H 2 L and β N H 2 are set to be 0.0025 and 0.261, respectively, which are similar to those in [36] for the Virginia opossum. The third type of host species is non-competent, β H 2 L = β N H 2 =0, such as the western fence lizard. The authors in [44, 45] stated that the western fence lizard is not able to spread the Lyme-pathogen since the species has a powerful immune system so that it can clean up the Lyme-pathogen when it is bitten by an infected tick. The death rate of each host species is set as μ H 2 = 0.0027 per day due to their similar life spans.

Time-dependent parameters

In order to investigate the impact of climate warming on the seasonal tick population abundance and Lyme-pathogen invasion, temperature is considered as a variable index in our study and we assume that it changes periodically with time. Therefore, those time-dependent coefficients are indeed temperature-dependent, and periodic in time. In order to parameterize these coefficients, we first estimate these values at a discrete manner at each day of a year, then these coefficients are smoothed into a continuous manner by employing Fourier series. In the remaining of this subsection, we will estimate each time-dependent coefficient at each day of a year.

To begin with this, the model is parameterized for the location Long Point, reported to be the first tick endemic area in Canada. Two temperature datasets for this area are collected from nearby meteorological stations, the Port Dover for the period 1961−1990 and Delhi CDA for the period 1981−2010 due to the unavailability of Port Dover Station recently. For these two stations, the 30-year normal temperature data are collected from the Environment Canada website (Figure 2) [46].
Figure 2

30 year normal mean monthly temperature under two settings near Long Point. The blue solid and red dashed curves represent the monthly temperature for the periods 1961−1990 period and 1981−2010, respectively. We set monthly temperature to be 0°C if it is lower than 0°C. Both are collected from Environment Canada website [46].

Next we turn to the estimation of time-dependent development rates: b(t i ), d E (t i ), d L (t i ) and d N (t i ), at day t i for t i =1,2,,365. To estimate these values, the following relations [1, 8, 4749]
D 1 ( t i ) = 1300 × T ( t i ) 1.42 + 10 ( pre-oviposition period ) ,
(6)
D 2 ( t i ) = 34234 × T ( t i ) 2.27 + 21 ( egg to larva period ) ,
(7)
D 3 ( t i ) = 101181 × T ( t i ) 2.55 + 3 ( larva to nymph period ) ,
(8)
D 4 ( t i ) = 1596 × T ( t i ) 1.21 + 5 ( nymph to adult period )
(9)

will be used, where T(t i ) represents temperature at the specific day t i in unit Celsius (°C). Using the same method presented in [2], the birth rate b(t i ) is directly obtained from the product of maximum number of eggs p produced and the reciprocal of duration of pre-oviposition period at day t i as shown Eq. 6, namely b(t i )=p/D1(t i ). The development rate of eggs d E (t i ) is directly calculated as reciprocal of development duration from egg to larva (Eq.7). The calculation of development rate of nymphs d N (t i ) is composed by two cases: (i) it is directly estimated as a reciprocal of development duration from nymph to adult (Eq. 9) before diapause; (ii) it is calculated by the method in [2] during diapause. The estimate of larval development rate d L (t i ) is a bit complex. We first consider the concept of the daily development proportion of larvae which is calculated as the reciprocal of development duration from larva to nymph at some specific days (Eq. 8) [2]. To obtain d L (t i ), we therefore calculate all daily development proportions from day t i until day t i +n for the subsequent n days such that the sum of these proportions reaches unity, then n is regarded as the development duration of larvae at the specific day t i . Finally, d L (t i ) is estimated as 1 n which is dependent of the temperatures of subsequent days.

The feeding rates F L (t i ), F N (t i ) and F A (t i ), affected by both hosts abundance and ambient temperatures, are directly calculated from the following formulas [1]:
F L ( t i ) = 0.0013 ( H 1 + p 1 H 2 ) 0.515 θ L ( T ( t i ) ) , F N ( t i ) = 0.0013 ( H 1 + p 2 H 2 ) 0.515 θ N ( T ( t i ) ) , F A ( t i ) = 0.086 D 0.515 θ A ( T ( t i ) ) ,
where θ L (T(t i )), θ N (T(t i )) and θ A (T(t i )) represent questing activity proportions at respective tick stage at day t i which are calibrated with data from Public Health Agency of Canada (personal communication). We refer the readers to the literature [2] for more details on the estimation of these periodic parameters. Figure 3 shows the patterns of these time-dependent parameters in one-year for the case p1= p2=0.
Figure 3

Development rates and feeding rates of I. scapularis ticks within one year period. The blue solid and red dashed curves are related to the associated development rates and feeding rates under temperatures in the periods 1961−1990 and 1981−2010, respectively; The numbers at the left top corner in each subfigure indicate the areas under the associated curves, which are used to differentiate the differences of these rates under the two temperature settings.

Results

We use various indices to measure the Lyme disease risk to humans: (i) R v , used to determine the tick population persistence; (ii) R 0 , as an index for the pathogen population persistence; (iii) density of questing nymphs (DON) in a seasonal pattern; (iv) density of infected questing nymphs (DIN), which reveals the absolute risk of Lyme disease by showing the absolute amount of infected ticks and the pattern of seasonality; and (v) nymphal infection prevalence (NIP) in a seasonal pattern, the proportion of the number of infected questing nymphs in total number of questing nymphs, which characterizes the degree of humans to be infected. All these are widely used indices and we use them to jointly measure the Lyme disease risk to humans [1, 12, 16, 18, 26, 28, 31].

In all simulations, every solution, irrespective of the initial values, of the model system (2) approaches to a seasonal state which is consistent with the theoretical results. Moreover, disease risk goes extinct when R 0 < 1 , while the seasonal risk pattern appears when R 0 > 1 . The numerical calculation of R v is implemented by the dichotomy method where the system d X/d t= ( F ( t ) / R v V ( t ) ) X ( t ) has a dominant Floquet multiplier equal to 1 [50]. A similar method is used to estimate R 0 . In what follows, all results are based on the model outputs at the steady state by running 40 years simulations.

Impact of climate warming on tick population growth and pathogen transmission

To study the potential effect of climate warming on disease risk, we compare simulations for two different temperature settings, at periods 1961−1990 and 1981−2010, with the absence of alternative host species. The curves of time-dependent parameters under these two temperature settings are shown in Figure 3. Moreover, the numbers on the upper left corner represent the areas under the corresponding curves, reflecting the variation of time-dependent parameters in different temperature consitions. We notice that the development rates and the feeding rates of immature ticks increase with increased temperature. However the feeding rate of adults decreases instead, which is because adult ticks have the limiting host seeking capacity when the temperature is too low or high [1].

With climate warms up from the period 1961−1990 to 1981−2010, the value of R v increases from 1.38 to 1.62, and the values of R 0 also increases from 0.90 (below unity) to 1.19 (above unity). As shown in Figure 4, our simulations confirm the persistence of tick population when R v > 1 and establishment of pathogen population if R 0 > 1 . These are in agreement with the theoretical conclusions. We also notice that the number of questing nymphs increases with higher temperature (Figures 4(a), (c)). Moreover, the pattern of infected questing nymphs changes from extinction to an absolutely positive stable oscillation showing the emergence of disease risk (Figures 4(b), (d)). It is important to notice that the active window of (infected) nymphs has been slightly enlarged with warmer temperature (Figure 4(c)). In summary, our study shows that climate warming plays an important role to accelerate the reproduction of the tick population and extend their active windows, and therefore increase the risk of Lyme disease. Moreover, the pattern of seasonality for ticks and pathogens may be changed with the temperature.
Figure 4

The variations in the sizes of total questing nymphs and infected questing nymphs with the two temperature settings mentioned above. The red solid curves represent the outputs by seeding the model with 1961−1990 temperature data ( R v = 1.38 and R 0 = 0.90 in this case), while the blue dashed curves represent the model outputs by 1981−2010 temperature data ( R v = 1.62 and R 0 = 1.19 in this case). (a) Total questing nymphs; (b) infected questing nymphs in the 40 year simulation; (c) seasonality of questing nymphs at the steady state; (d) seasonality of infected questing nymphs at the steady state, where shaded portions in both (c) and (d) represent the active seasons of the questing nymphs.

Impact of host biodiversity on the disease risk

Now, we seed the model with temperature condition in the 1981−2010 period so that time-dependent birth rate and development rates remain the same. However, we add an alternative host species to the original host community which is assumed to be composed of the white-footed mice and the white-tailed deer alone. This permits us to study the potential impact of host biodiversity on the risk of Lyme disease. Then, the number of the alternative host species will change the density-dependent death rates and the feeding rates of ticks.

As shown in Figure 5, regardless of the newly introduced alternative species, we always observe that the values of R v continuously increase with the increased number of hosts; while the change of R 0 is closely connected to the species of the introduced hosts. Introduction of new hosts will always provide more food for the ticks and thus promotes the growth of tick population. However, the variation of the disease risk is not as simple as we imagine. For instance, the values of R 0 persistently increase with the increased number of the eastern chipmunk introduced, however continuously decrease for the Virginia opossum, while first increase then decrease for the western fence lizard (Figure 5). For the eastern chipmunk, recognized as the type with a high reservoir competence ( β H 2 L = 0.569 and β N H 2 = 0.971 ), their ability of Lyme-pathogen transmission and high biting bias coefficient of nymphs (p2=3.5) facilitate the growth of tick population and spread of the pathogen. For the Virginia opossum with a low reservoir competence ( β H 2 L = 0.0025 and β N H 2 = 0.261 ), the reduction of R 0 largely attributes to not only the low transmission ability, but also their large biting biases coefficients (p1=7.2, p2=36.9). In this scenario, a great amount of tick bites are attracted to the low competent hosts, and infectious bites are wasted on this incompetent host. For the case of the western fence lizard, we also observe that R 0 increases at the small size of this species even it is a non-competent host, but eventually reduces when the size of western fence lizard attains a certain level.
Figure 5

Log plots of variations of ratios R v and R 0 against the number of alternative hosts. In the case of the eastern chipmunks, p1=0.4, p2=3.5, β H 2 L = 0.569 , β N H 2 = 0.971 , μ H 2 = 0.0027 ; for the western fence lizard, p1=1, p2=1, β H 2 L = 0 , β N H 2 = 0 , μ H 2 = 0.0027 ; for the Virginia opossum, p1=7.2, p2=36.9, β H 2 L = 0.0025 , β N H 2 = 0.261 , μ H 2 = 0.0027 . For all simulations, the temperature condition is fixed on the period 1981−2010.

To clearly understand the “dilute effect” and “amplification effect” in this respect, we would like to examine three indices: DON, DIN and NIP. As shown in Figure 6, introduction of numbers of the eastern chipmunk from 10, 20 to 40 leads to continuous increase of DON, DIN and NIP, and this indicates that the eastern chipmunk offers an efficient host species to amplify the risk of Lyme disease; if the same numbers of the Virginia opossum as these of the eastern chipmunk are added, we notice that DON increases, but both DIN and NIP decrease instead, input of this species indeed reflects the “dilute effect” through reducing not only the absolute amount of infected ticks, but also the proportion of infection; we are surprised to observe that DON continuously increases, DIN first increases and then decreases, while NIP continuously decreases when the western fence lizard is added into the existing host community. That is, this non-competent additional species amplifies the risk of Lyme disease in the sense of absolute amount; on the contrary it also dilutes the risk in the sense of relative proportion of infection. This finding is in good agreement with the debate raised in [14, 51], where authors revealed that the western fence lizard, as a non-competent host, does not always dilute the risk of Lyme disease.
Figure 6

Variations in the sizes of DON, DIN and NIP under different host sizes. The numbers on the left panel indicate that the sizes of associated alternative hosts are added into the host community. The scenarios where the eastern chipmunk is added are shown on the upper panel. The middle panel shows the scenarios where the Virginia opossum is considered as the alternative host; the bottom panel shows the situations where the western fence lizard is added. All the associated parameter values are the same as those in Figure 5 except the sizes of alternative hosts.

We also perform sensitivity analysis of the threshold R 0 against the biting biases p1 and p2. The result shows that R 0 is very sensitive to the variations of both biting biases (Figure 7). Moreover, the relationships between R 0 and p1 and p2 varies with host species: (i) R 0 increases with increased p1 and p2 in the case of the eastern chipmunk, and therefore this species always facilitates disease transmission within our parameter region; (ii) the relation between R 0 and the larvae bias p1 is neither positive nor negative for the case of the western fence lizard or the Virginia opossum, which implies both “dilution effect” and “amplification effect” would occur.
Figure 7

Relationship between the reproduction ratio R 0 , p 1 and p 2 for the three types of alternative hosts: the eastern chipmunk, the Virginia opossum and the western fence lizard. The number of alternative hosts is set as 30 and all other parameter values are same as those in Figure 5.

Discussion and conclusion

In this paper, we developed a periodic deterministic system of ordinary differential equations to investigate the impact of both climate condition and host biodiversity on Lyme disease pathogen transmission through the mathematical analysis and computer simulations. The model was parameterized using field and local ecological and epidemiological data. The critical ratios, R v and R 0 , in combination with other widely used indices, can then provide pivotal information on the impact of temperature variation and host biodiversity on Lyme disease spread.

We found that climate warming facilitates the reproduction of I. scapularis population and accelerates the spread of Lyme-pathogen, and then increases the risk of Lyme disease infection. Furthermore, we also have noticed that climate change can slightly change the seasonality of the infected questing nymphs and slightly broaden the active period of the infected questing nymphs, and therefore slightly change the seasonality of the risk of Lyme disease. However, when a new host species was added, we didn’t observe the change of seasonality of the tick population, but we observed the increase of the quantity of total ticks including infected ticks.

The impact of host biodiversity on the Lyme disease risk is a complex issue and remains challenging in conservation ecology and zoonotic epidemiology. However, this issue has both theoretical and practical importance since this may reveal whether the biodiversity conservation can be used as an effective measure for the prevention and control of the zoonotic disease. For Lyme disease, both the dilution effect [5, 5257] and amplification effect [14] have been observed through field and theoretical studies, where many factors such as spatial scale, host competition, host resistance, tick contact rate were considered [26, 37, 58, 59]. Through this modeling study, both “amplification effect” and “dilution effect” have been observed, where multiple indices ( R v , R 0 , DON, DIN and NIP) instead of a single index were utilized. However, the effect does not depend upon the host competence alone, but is a joint outcome of current climate condition, host transmission ability, the numbers of hosts and so on.

In conclusion, climate warming plays a crucial role to speed up the spread of Lyme disease and hence increase the disease risk since climate warming can promote the tick population growth. Introduction of new host species into host community can certainly increase the amount of total ticks, but is not necessary increase the number of infected ticks. In order to obtain a definitive answer to the question “How does the biodiversity of the host community affect the disease risk?”, reliable field study in combination with local abiotic and biotic factors is necessary.

By assuming a spatially homogeneous habitat, the model formulated here has not evaluated the effect of spatial heterogeneity on disease pattern. As ticks can disperse mainly due to its host movement, such as short distance movement due to rodents, long distance travel due to deer [18] and even longer distance because of the bird migration [60]. In 2002, Caraco et al. [18] proposed a reaction-diffusion model for Lyme disease in the northeast United States to investigate the spreading speed of the Lyme disease. The global dynamics of this model was further anlyzed in [61]. A periodic reaction-diffusion system was proposed to study the impact of spatial structure and seasonality on the spreading of the pathogen [31]. The effect of bird migration on Lyme dispersal was studied in [62]. It would be interesting to incorporate our current model formulation into the aforementioned studies involving spatial aspect of Lyme disease spread to address the complicated spatiotemporal spread patterns of Lyme disease with biodiversity and seasonal variation.

Appendix 1: Proof of Theorem 2.

Proof.

It follows from [63], Theorem 5.2.1, that for any initial value x0X, system (1) admits a unique nonnegative solution x(t,x0) through this initial value with the maximal interval of existence [0,σ) for some σ>0.

Let L F =L FS +L FI , N Q =N QS +N QI , N F =N FS +N FI , A Q =A QS +A QI and A F =A FS +A FI . Then we can see that the tick growth is governed by the following system:
dE dt = b ( t ) A F ( t ) μ E E ( t ) d E ( t ) E ( t ) , d L Q dt = d E ( t ) E ( t ) μ QL L Q ( t ) F L ( t ) L Q ( t ) , d L F dt = F L ( t ) L Q ( t ) μ FL L F ( t ) D L L F 2 ( t ) d L ( t ) L F ( t ) , d N Q dt = d L ( t ) L F ( t ) μ QN N Q ( t ) F N ( t ) N Q ( t ) , d N F dt = F N ( t ) N Q ( t ) μ FN N F ( t ) D N N F 2 ( t ) d N ( t ) N F ( t ) , d A Q dt = d N ( t ) N F ( t ) μ QA A Q ( t ) F A ( t ) A Q ( t ) , d A F dt = F A ( t ) A Q ( t ) μ FA A F ( t ) D A A F 2 ( t ) .
(10)
For any periodic nonnegative function f(t) with period τ, denote f ̂ = max t [ 0 , τ ] f ( t ) and f ~ = min t [ 0 , τ ] f ( t ) . It is easy to see that system (10) can be controlled by the following cooperative system:
d u 1 dt = b ̂ u 7 ( t ) μ E u 1 ( t ) , d u 2 dt = d E ̂ u 1 ( t ) μ QL u 2 ( t ) , d u 3 dt = F L ̂ u 2 ( t ) μ FL u 3 ( t ) , d u 4 dt = d L ̂ u 3 ( t ) μ QN u 4 ( t ) , d u 5 dt = F N ̂ u 4 ( t ) μ FN u 5 ( t ) , d u 6 dt = d N ̂ u 5 ( t ) μ QA u 6 ( t ) , d u 7 dt = F A ̂ u 6 ( t ) μ FA u 7 ( t ) D A u 7 2 ( t ) .
(11)
Clearly, there is only one nonnegative equilibrium zero for system (11) when
F A ̂ d N ̂ μ QA F N ̂ μ FN d L ̂ μ QN F L ̂ μ FL d E ̂ μ QL b ̂ μ E μ FA .

If F A ̂ d N ̂ μ QA F N ̂ μ FN d L ̂ μ QN F L ̂ μ FL d E ̂ μ QL b ̂ μ E > μ FA , system (11) admits another positive equilibrium. It then follows from [64], Corollary 3.2, that either zero is globally asymptotically stable or the positive equilibrium is globally asymptotically stable for all nonzero solutions. Hence the comparison principle implies that (E(t), L Q (t), L F (t), N Q (t), N F (t), A Q (t), A F (t)) is bounded for any t[0,σ). Thus, we see that σ= and the solution for model (1) is bounded and exists globally for any nonnegative initial value. □

Appendix 2: Proof of Theorem 2.2

Proof.

Theorem 2.3.4 in [42] directly implies that if R v 1 , then zero is globally asymptotically stable for system (3) in + 7 ; if R v > 1 , then system (3) admits a unique 6τ-positive periodic solution
( E ( t ) , L Q ( t ) , L F ( t ) , N Q ( t ) , N F ( t ) , A Q ( t ) , A F ( t ) ) ,
and it is globally asymptotically stable for system (3) with initial values in + 7 { 0 } . It remains to prove that the 6τ-positive periodic solution (E(t), L Q ( t ) , L F ( t ) , N Q ( t ) , N F ( t ) , A Q ( t ) , A F ( t ) ) is also τ-periodic. Since for any x + 7 { 0 } , lim n P 6 n ( x ) = (E(0), L Q ( 0 ) , L F ( 0 ) , N Q ( 0 ) , N F ( 0 ) , A Q ( 0 ) , A F ( 0 ) ) where P is the Poincar e ́ map associated with the τ-periodic system (3). Hence,
lim n P ( P 6 n ( x ) ) = P ( E ( 0 ) , L Q ( 0 ) , L F ( 0 ) , N Q ( 0 ) , N F ( 0 ) , A Q ( 0 ) , A F ( 0 ) ) .
On the other hand,
lim n P 6 n ( P ( x ) ) = ( E ( 0 ) , L Q ( 0 ) , L F ( 0 ) , N Q ( 0 ) , N F ( 0 ) , A Q ( 0 ) , A F ( 0 ) ) .
Thus,
P ( E ( 0 ) , L Q ( 0 ) , L F ( 0 ) , N Q ( 0 ) , N F ( 0 ) , A Q ( 0 ) , A F ( 0 ) ) = ( E ( 0 ) , L Q ( 0 ) , L F ( 0 ) , N Q ( 0 ) , N F ( 0 ) , A Q ( 0 ) , A F ( 0 ) ) ,

which implies that (E(t), L Q ( t ) , L F ( t ) , N Q ( t ) , N F ( t ) , A Q ( t ) , A F ( t ) ) is τ-periodic. □

Appendix 3: Proof of Theorem 2.

Proof.

We first consider the τ-periodic system as a 11τ-periodic system. Let P be the Poincar e ́ map of system (2), that is, P(x0)=x(11τ,x0), where x(t,x0) is the solution of system (2) through x0. Then P is compact. Let ω = ω(x0) be the omega limit set of P(x0). It then follows from [65], Lemma 2.1, (see also [42], Lemma 1.2.1) that ω is an internally chain transitive set for P.

(i) In the case where R v 1 , we obtain lim t x i ( t ) = 0 for i[1,9]. Hence, ω= {(0, 0, 0, 0, 0, 0, 0, 0, 0)}×ω 1 for some ω 1 2 . It is easy to see that
P | ω ( 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , H 1 I ( 0 ) , H 2 I ( 0 ) ) = ( 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , P 1 ( H 1 I ( 0 ) , H 2 I ( 0 ) ) ) ,
where P1 is the Poincar e ́ map associated with the following equation:
d H 1 I dt = μ H 1 H 1 I , d H 2 I dt = μ H 2 H 2 I .
(12)

Since ω is an internally chain transitive set for P, it easily follows that ω1 is an internally chain transitive set for P1. Since {0} is globally asymptotically stable for system (12), [65], Theorem 3.2, implies that ω1={(0,0)}. Thus, we have ω={0}, which proves that every solution converges to zero.

(ii) In the case where R v > 1 , then there exists a positive periodic solution, (E (t), L Q ( t ) , L F ( t ) , N Q ( t ) , N F ( t ) , A Q ( t ) , A F ( t ) ), for system (3) such that for any x 0 with i = 1 7 x i 0 > 0 , we have
lim t { ( x 1 ( t ) , x 2 ( t ) , x 3 ( t ) , x 4 ( t ) , x 5 ( t ) , x 6 ( t ) , x 7 ( t ) ) ( E ( t ) , L Q ( t ) , L F ( t ) , N Q ( t ) , N F ( t ) , A Q ( t ) , A F ( t ) ) } = 0 .
Thus, ω = { ( E ( 0 ) , L Q ( 0 ) , L F ( 0 ) , N Q ( 0 ) , N F ( 0 ) , A Q ( 0 ) , A F ( 0 ) ) } × ω 2 for some ω 2 4 , and
P | ω E ( 0 ) , L Q ( 0 ) , L F ( 0 ) , N Q ( 0 ) , N F ( 0 ) , A Q ( 0 ) , A F ( 0 ) , x 8 , x 9 , x 10 , x 11 = E ( 0 ) , L Q ( 0 ) , L F ( 0 ) , N Q ( 0 ) , N F ( 0 ) , A Q ( 0 ) , A F ( 0 ) , P 2 x 8 , x 9 , x 10 , x 11 ,
where P2 is the Poincaré map associated with system (5). Since ω is an internally chain transitive set for P, ω2 is an internally chain transitive set for P2. Since R 0 1 , {(0,0,0,0)} is globally asymptotically stable for system (5) according to Theorem 2.3. It then follows from [65], Theorem 3.2, that ω2={0}. This proves
ω = { ( E ( 0 ) , L Q ( 0 ) , L F ( 0 ) , N Q ( 0 ) , N F ( 0 ) , A Q ( 0 ) , A F ( 0 ) , 0 , 0 , 0 , 0 ) } .

Therefore, statement (ii) holds.

(iii) In the case where R v > 1 and R 0 > 1 , then there exists a positive periodic solution, ( E ( t ) , L Q ( t ) , L F ( t ) , N Q ( t ) , N F ( t ) , A Q ( t ) , A F ( t ) ) , for system (3) such that for any x 0 with i = 1 7 x i 0 > 0 , we have
lim t ( x 1 ( t ) , x 2 ( t ) , x 3 ( t ) , x 4 ( t ) , x 5 ( t ) , x 6 ( t ) , x 7 ( t ) ) ( E ( t ) , L Q ( t ) , L F ( t ) , N Q ( t ) , N F ( t ) , A Q ( t ) , A F ( t ) ) = 0 .
It then follows that ω = { ( E ( 0 ) , L Q ( 0 ) , L F ( 0 ) , N Q ( 0 ) , N F ( 0 ) , A Q ( 0 ) , A F ( 0 ) ) ) } × ω 3 for some ω 3 4 , and
P | ω E ( 0 ) , L Q ( 0 ) , L F ( 0 ) , N Q ( 0 ) , N F ( 0 ) , A Q ( 0 ) , A F ( 0 ) , x 8 , x 9 , x 10 , x 11 = E ( 0 ) , L Q ( 0 ) , L F ( 0 ) , N Q ( 0 ) , N F ( 0 ) , A Q ( 0 ) , A F ( 0 ) , P 2 x 8 , x 9 , x 10 , x 11 ,

where P2 is the solution semiflow of system (5). Since ω is an internally chain transitive set for P, it follows that ω3 is an internally chain transitive set for P2. We claim that ω3≠{0} for any ( x 8 0 , x 9 0 , x 10 0 , x 11 0 ) > 0 .

Assume that, by contradiction, ω3={0}. That is
ω = { ( E ( 0 ) , L Q ( 0 ) , L F ( 0 ) , N Q ( 0 ) , N F ( 0 ) , A Q ( 0 ) , A F ( 0 ) , 0 , 0 , 0 , 0 ) }
for some ( x 8 0 , x 9 0 , x 10 0 , x 11 0 ) > 0 . Then, we have
lim t ( x ( t ) ( E ( t ) , L Q ( t ) , L F ( t ) , N Q ( t ) , N F ( t ) , A Q ( t ) , A F ( t ) , 0 , 0 , 0 , 0 ) ) = 0 .
(13)
Since R 0 > 1 , there exists some δ>0 such that the spectral radius of the Poincar e ́ map associated with the following linearized system is greater than unity:
d L FI dt = ( β H 1 L H 1 I ( t ) H 1 + p 1 H 2 + β H 2 L p 1 H 2 I ( t ) H 1 + p 1 H 2 ) F L ( t ) ( L Q ( t ) δ ) D L ( L F ( t ) + δ ) L FI ( t ) ( d L ( t ) + μ FL ) L FI ( t ) , d N QI dt = d L ( t ) L FI ( t ) ( μ QN + F N ( t ) ) N QI ( t ) , d H 1 I dt = F N ( t ) β N H 1 N QI ( t ) H 1 H 1 I ( t ) H 1 + p 2 H 2 μ H 1 H 1 I ( t ) , d H 2 I dt = F N ( t ) β N H 2 N QI ( t ) p 2 ( H 2 H 2 I ( t ) ) H 1 + p 2 H 2 μ H 2 H 2 I ( t ) .
It then follows from the same argument as in the proof of Theorem 2.3 that the following system
d u 1 dt = β H 1 L u 3 ( t ) H 1 + p 1 H 2 + β H 2 L p 1 u 4 ( t ) H 1 + p 1 H 2 F L ( t ) L Q ( t ) δ D L L F ( t ) + δ u 1 ( t ) d L ( t ) + μ FL u 1 ( t ) , d u 2 dt = d L ( t ) u 1 ( t ) μ QN + F N ( t ) u 2 ( t ) , d u 3 dt = F N ( t ) β N H 1 u 2 ( t ) H 1 u 3 ( t ) H 1 + p 2 H 2 μ H 1 u 3 ( t ) , d u 4 dt = F N ( t ) β N H 2 u 2 ( t ) p 2 H 2 u 4 ( t ) H 1 + p 2 H 2 μ H 2 u 4 ( t )
admits a positive periodic u(t) such that
lim t ( u ( t ) u ( t ) ) = 0 , u ( t ) 0 .
Therefore, there exists some τ0>0 such that for all t>τ0,
( x 1 ( t ) , , x 7 ( t ) ) ( E ( t ) , L Q ( t ) , L F ( t ) , N Q ( t ) , N F ( t ) , A Q ( t ) , A F ( t ) ) δ.
Hence, we conclude that
d x 8 dt β H 1 L x 10 ( t ) H 1 + p 1 H 2 + β H 2 L p 1 x 11 ( t ) H 1 + p 1 H 2 F L ( t ) L Q ( t ) δ D L L F ( t ) + δ x 8 ( t ) d L ( t ) + μ FL x 8 ( t ) , d x 9 dt = d L ( t ) x 8 ( t ) μ QN + F N ( t ) x 9 ( t ) , d x 10 dt = F N ( t ) β N H 1 x 9 ( t ) H 1 x 10 ( t ) H 1 + p 2 H 2 μ H 1 x 10 ( t ) , d x 11 dt = F N ( t ) β N H 2 x 9 ( t ) p 2 H 2 x 11 ( t ) H 1 + p 2 H 2 μ H 2 x 11 ( t ) ,
for all t>τ0. By a standard comparison argument, we have
liminf t x 8 ( t ) , x 9 ( t ) , x 10 ( t ) , x 11 ( t ) u ( t ) 0 ,

a contradiction to (13).

Since ω3≠{0} and the positive periodic solution ( L FI ( t ) , N QI ( t ) , H 1 I ( t ) , H 2 I ( t ) ) is globally asymptotically stable for system (5) in + 4 { 0 } , it follows that
ω 3 W s ( ( L FI ( 0 ) , N QI ( 0 ) , H 1 I ( 0 ) , H 2 I ( 0 ) ) ) ,
where W s L FI ( 0 ) , N QI ( 0 ) , H 1 I ( 0 ) , H 2 I ( 0 ) is the stable set for ( L FI ( 0 ) , N QI ( 0 ) , H 1 I ( 0 ) , H 2 I ( 0 ) ) with respect to the Poincar e ́ map P2. By [65], Theorem 3.1, we then get
ω 3 = L FI ( 0 ) , N QI ( 0 ) , H 1 I ( 0 ) , H 2 I ( 0 ) .
Thus,
ω = E ( 0 ) , L Q ( 0 ) , L F ( 0 ) , N Q ( 0 ) , N F ( 0 ) , A Q ( 0 ) , A F ( 0 ) , L FI ( 0 ) , N QI ( 0 ) , H 1 I ( 0 ) , H 2 I ( 0 ) ,

and hence, statement (iii) is valid.

At last, using a similar argument as in the proof of Theorem 2.2, we can show that the globally attractive 11τ-periodic solution in each case is also τ-periodic solution. □

Declarations

Acknowledgements

This work was partially supported by the Mitacs, the GEOIDE, NSFC (11301442) and RGC (PolyU 253004/14P).

Authors’ Affiliations

(1)
Department of Applied Mathematics, The Hong Kong Polytechnic University
(2)
Department of Mathematics and Statistics, York University
(3)
Centre for Disease Modelling, York University
(4)
Department of Applied Mathematics, Western University

References

  1. Ogden NH, Bigras-Poulin M, O’Callaghan CJ, Barker IK, Lindsay LR, Maarouf A, Smoyer-Tomic KE, Waltner-Toews D, Charron D: A dynamic population model to investigate effects of climate on geographic range and seasonality of the tick Ixodes scapularis. Int J Parasitol. 2005, 35: 375-389. 10.1016/j.ijpara.2004.12.013.View ArticlePubMedGoogle Scholar
  2. Wu X, Duvvuri VR, Lou Y, Ogden NH, Pelcat Y, Wu J: Developing a temperature-driven map of the basic reproductive number of the emerging tick vector of Lyme disease Ixodes scapularis in Canada. J Theor Biol. 2013, 319: 50-61.View ArticlePubMedGoogle Scholar
  3. CDC: Summary of notifiable diseases–United States, 2010. Morb Mortal Wkly Rep. 2012, 59: 1-111.Google Scholar
  4. Ostfeld RS: Lyme Disease: The Ecology of a Complex System. 2011, New York: Oxford University PressGoogle Scholar
  5. LoGiudice K, Ostfeld RS, Schmidt KA, Keesing F: The ecology of infectious disease: effects of host diversity and community composition on Lyme disease risk. Proc Natl Acad Sci USA. 2003, 100: 567-571. 10.1073/pnas.0233733100.PubMed CentralView ArticlePubMedGoogle Scholar
  6. Duffy DC, Campbell SR: Ambient air temperature as a predictor of activity of adult Ixodes scapularis (Acari: Ixodidae). J Med Entomol. 1994, 31: 178-180. 10.1093/jmedent/31.1.178.View ArticlePubMedGoogle Scholar
  7. Loye JE, Lane RS: Questing behavior of Ixodes pacificus (Acari: Ixodidae) in relation to meteorological and seasonal factors. J Med Entomol. 1988, 25: 391-398. 10.1093/jmedent/25.5.391.View ArticlePubMedGoogle Scholar
  8. Ogden NH, Lindsay LR, Beauchamp G, Charron D, Maarouf A, O’Callaghan CJ, Waltner-Toews D, Barker IK: Investigation of relationships between temperature and developmental rates of tick Ixodes scapularis (Acari: Ixodidae) in the laboratory and field. J Med Entomol. 2004, 41: 622-633. 10.1603/0022-2585-41.4.622.View ArticlePubMedGoogle Scholar
  9. Vail SC, Smith G: Vertical movement and posture of blacklegged tick (Acari: Ixodidae) nymphs as a function of temperature and relative humidity in laboratory experiments. J Med Entomol. 2002, 39: 842-846. 10.1603/0022-2585-39.6.842.View ArticlePubMedGoogle Scholar
  10. Ghosh M, Pugliese A: Seasonal population dynamics of ticks, and its influence on infection transmission: a semi-discrete approach. Bull Math Biol. 2004, 66: 1659-1684. 10.1016/j.bulm.2004.03.007.View ArticlePubMedGoogle Scholar
  11. Ogden NH, Radojević M, Wu X, Duvvuri VR, Leighton P, Wu J: Estimated effects of projected climate change on the basic reproductive number of the Lyme disease vector Ixods scapularis. Environ Health Perspect. 2014, 122: 631-638.PubMed CentralPubMedGoogle Scholar
  12. Norman R, Bowers RG, Begon M, Hudson PJ: Persistence of tick-borne virus in the presence of multiple host species: tick reservoirs and parasite-mediated competition. J Theor Biol. 1999, 200: 111-118. 10.1006/jtbi.1999.0982.View ArticlePubMedGoogle Scholar
  13. Ogden NH, Bigras-Poulin M, O’Callaghan CJ, Barker IK, Kurtenbach K, Lindsay LR, Charron DF: Vector seasonality, host infection dynamics and fitness of pathogens transmitted by the tick Ixodes scapularis. Parasitology. 2007, 134: 209-227. 10.1017/S0031182006001417.View ArticlePubMedGoogle Scholar
  14. Randolph SE, Dobson ADM: Pangloss revisited: a critique of the dilution effect and the biodiversity-buffers-disease paradigm. Parasitology. 2012, 139: 847-863. 10.1017/S0031182012000200.View ArticlePubMedGoogle Scholar
  15. Rosá R, Pugliese A: Effects of tick population dynamics and host densities on the persistence of tick-borne infections. Math Biosci. 2007, 208: 216-240. 10.1016/j.mbs.2006.10.002.View ArticlePubMedGoogle Scholar
  16. Rosá R, Pugliese A, Norman R, Hudson PJ: Thresholds for disease persistence in models for tick-borne infections including non-viraemic transmission, extended feeding and tick aggregation. J Theor Biol. 2003, 224: 359-376. 10.1016/S0022-5193(03)00173-5.View Article