# Impact of biodiversity and seasonality on Lyme-pathogen transmission

- Yijun Lou
^{1}, - Jianhong Wu
^{2, 3}Email author and - Xiaotian Wu
^{3, 4}

**11**:50

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

© Lou et al.; licensee BioMed Central Ltd. 2014

**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

## 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 [6–9], 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 [12–16], stage structure of ticks [1, 13, 17–22] 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.

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

*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.

**Variable explanations used in the model (1)**

Variable | Meaning |
---|---|

| the number of eggs |

| the number of questing larvae |

| the number of susceptible feeding larvae |

| the number of infected feeding larvae |

| the number of susceptible questing nymphs |

| the number of infected questing nymphs |

| the number of susceptible feeding nymphs |

| the number of infected feeding nymphs |

| the number of susceptible questing adults |

| the number of infected questing adults |

| the number of susceptible feeding adults |

| the number of infected feeding adults |

| the number of infected white-footed mice |

| 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 *H*_{1} (mainly *Peromyscus leucopus*) with the mortality rate ${\mu}_{{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 *H*_{2} with mortality rate ${\mu}_{{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 *H*_{1I} and *H*_{2I}, respectively. The Lyme-pathogen transmission cycle between the hosts and multi-stage tick population is presented in the diagram of Figure 1.

*p*

_{1}(

*p*

_{2}), to describe larval (nymphal) ticks biting bias on their hosts [33, 34]. Specifically,

*p*

_{1}>1(

*p*

_{2}>1) indicates one host

*H*

_{2}can attract more larval (nymphal) bites than one host

*H*

_{1}and vice versa when 0<

*p*

_{1}<1(0<

*p*

_{2}<1). Using the method described in [35], ${F}_{L}(t)\frac{{H}_{1}}{{H}_{1}+{p}_{1}{H}_{2}}\frac{{H}_{1I}(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 ${\beta}_{{H}_{1}L}{F}_{L}(t)\frac{{H}_{1}}{{H}_{1}+{p}_{1}{H}_{2}}\frac{{H}_{1I}(t)}{{H}_{1}}$ is the average infection rate at which a susceptible larva gets infected from mice, where ${\beta}_{{H}_{1}L}$ is the pathogen transmission probability per bite from infectious mice

*H*

_{1}to susceptible larvae. Using the same idea, the infection rate of larvae from the infected alternative host

*H*

_{2}can be accounted. Therefore, the larval infection rate is given by

*τ*-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.

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

Parameter | Meaning | (Value, [reference]) or estimation |
---|---|---|

| mortality rate of eggs | (0.0025, [1]) |

| mortality rate of questing larvae | (0.006, [1]) |

| mortality rate of questing nymphs | (0.006, [1]) |

| mortality rate of questing adults | (0.006, [1]) |

| natural mortality rate of feeding larvae | (0.038, A) |

| natural mortality rate of feeding nymphs | (0.028, A) |

| natural mortality rate of feeding adults | (0.018, A) |

| the number of white-footed mice | (200, [1]) |

${\beta}_{{H}_{1}L}$ | transmission probability from | (0.6, [13]) |

${\beta}_{N{H}_{1}}$ | transmission probability from nymphs to | (1, [13]) |

${\mu}_{{H}_{1}}$ | death rate of the white-footed mice | (0.012, [13]) |

| the number of alternative host | (variable) |

${\beta}_{{H}_{2}L}$ | transmission probability from | (variable, [36]) |

${\beta}_{N{H}_{2}}$ | transmission probability from nymphs to | (variable, [36]) |

${\mu}_{{H}_{2}}$ | death rate of the alternative host | (variable) |

| the number of deer | (20, [1]) |

| the maximum number of eggs produced | (3000, [1]) |

| larval biting bias for host | (variable, [37]) |

| nymphal biting bias for host | (variable, [37]) |

| birth rate of eggs produced | (see subsection ‘Time-dependent parameters’) |

| development rate of eggs | (see subsection ‘Time-dependent parameters’) |

| development rate of larvae | (see subsection ‘Time-dependent parameters’) |

| development rate of nymphs | (see subsection ‘Time-dependent parameters’) |

| feeding rate of larvae | (see subsection ‘Time-dependent parameters’) |

| feeding rate of nymphs | (see subsection ‘Time-dependent parameters’) |

| feeding rate of adults | (see subsection ‘Time-dependent parameters’) |

| density-dependent mortality rate of feeding larvae | ($\frac{0.001084}{{H}_{1}+{p}_{1}{H}_{2}}$, E) |

| density-dependent mortality rate of feeding nymphs | ($\frac{0.001084}{{H}_{1}+{p}_{2}{H}_{2}}$, E) |

| density-dependent mortality rate of feeding adults | ($\frac{0.001084}{D}$, E) |

### 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**.

*x*(

*t*,

*x*

^{0}). Moreover, the solution

*x*(

*t*,

*x*

^{0})remains in

*X*for any

*t*≥0. Here, $x\in {\mathbb{R}}_{+}^{14}$ denotes a generic point with components

*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

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

*F*(

*t*)=(

*f*

_{ ij }(

*t*))

_{7×7}, where

*f*

_{1,7}(

*t*)=

*b*(

*t*) and

*f*

_{i,j}(

*t*)=0 if (

*i*,

*j*)≠(1,7), and

*V*(

*t*) =

*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*),

*t*≥

*s*, is the evolution operator of the linear periodic system $\frac{\mathit{\text{dy}}}{\mathit{\text{dt}}}=-V(t)y$. That is, for each $s\in \mathbb{R}$, the 7×7 matrix

*Y*(

*t*,

*s*) satisfies

*I*is the 7×7 identity matrix. Let

*C*

_{ τ }be the Banach space of all

*τ*-periodic functions from to ${\mathbb{R}}^{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

*t*≥

*s*. Hence,

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

*G*:

*C*

_{ τ }→

*C*

_{ τ }by

Then the spectral radius of *G* is defined as ${\mathcal{R}}_{v}:=\rho (G)$. In what follows, we call ${\mathcal{R}}_{v}$ as a threshold for tick population dynamics.

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

Note that the Poincar$\stackrel{\u0301}{\mathrm{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$\stackrel{\u0301}{\mathrm{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$\stackrel{\u0301}{\mathrm{e}}$ map associated with system (3) to obtain the following result, with the proof in Appendix 2.

**Theorem** **2.2**.

- (i)
If ${\mathcal{R}}_{v}\le 1$, then zero is globally asymptotically stable for system (3) in ${\mathbb{R}}_{+}^{7}$;

- (ii)If ${\mathcal{R}}_{v}>1$, then system (3) admits a unique
*τ*-positive periodic solution$({E}^{\ast}(t),{L}_{Q}^{\ast}(t),{L}_{F}^{\ast}(t),{N}_{Q}^{\ast}(t),{N}_{F}^{\ast}(t),{A}_{Q}^{\ast}(t),{A}_{F}^{\ast}(t)),$

and it is globally asymptotically stable for system (3) with initial values in ${\mathbb{R}}_{+}^{7}\setminus \left\{0\right\}$.

### The global dynamics of the full model

*t*,

*s*),

*t*≥

*s*, is the evolution operator of the linear periodic system $\frac{\mathit{\text{dy}}}{\mathit{\text{dt}}}=-\stackrel{~}{V}(t)y$. Let ${\stackrel{~}{C}}_{\tau}$ be the Banach space of all

*τ*-periodic functions from to ${\mathbb{R}}^{4}$, equipped with the maximum norm. Suppose $\varphi \in {\stackrel{~}{C}}_{\tau}$ is the initial distribution of infectious tick and host individuals in this periodic environment. Then $\stackrel{~}{F}(s)\varphi (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 $\stackrel{~}{Y}(t,s)\stackrel{~}{F}(s)\varphi (s)$ represents the distribution of those ticks who were newly produced at time

*s*and remain alive at time

*t*for

*t*≥

*s*. Hence,

*t*produced by all those infectious individuals

*ϕ*(

*s*) introduced at the previous time. Define the a next generation operator $\stackrel{~}{G}:{\stackrel{~}{C}}_{\tau}\to {\stackrel{~}{C}}_{\tau}$ by

It then follows from [39, 40] that the spectral radius of $\stackrel{~}{G}$ is define as ${\mathcal{R}}_{0}:=\rho (\stackrel{~}{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 ${\mathcal{R}}_{0}\le 1$, then zero is globally asymptotically stable for system (5) in ${\mathbb{R}}_{+}^{4}$; (ii) If ${\mathcal{R}}_{0}>1$, then system (5) admits a unique positive periodic solution $({L}_{\mathit{\text{FI}}}^{\ast}(t),{N}_{\mathit{\text{QI}}}^{\ast}(t),{H}_{1I}^{\ast}(t),{H}_{2I}^{\ast}(t))$ and it is globally asymptotically stable for system (5).

Based on the aforementioned two thresholds, ${\mathcal{R}}_{v}$ for ticks dynamics and ${\mathcal{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**.

*x*(

*t*,

*x*

^{0}) be the solution of system (2) through

*x*

^{0}. Then the following statements are valid:

- (i)
If ${\mathcal{R}}_{v}\le 1$, then zero is globally attractive for system (2);

- (ii)If ${\mathcal{R}}_{v}>1$ and ${\mathcal{R}}_{0}\le 1$, then$\begin{array}{l}\underset{\mathit{\text{t}}\to \infty}{lim}\{({x}_{1}(t),{x}_{2}(t),{x}_{3}(t),{x}_{4}(t),{x}_{5}(t),{x}_{6}(t),{x}_{7}(t))\\ \phantom{\rule{1em}{0ex}}-({E}^{\ast}(t),{L}_{Q}^{\ast}(t),{L}_{F}^{\ast}(t),{N}_{Q}^{\ast}(t),{N}_{F}^{\ast}(t),{A}_{Q}^{\ast}(t),{A}_{F}^{\ast}(t))\}=0,\end{array}$
and $\underset{\mathit{\text{t}}\to \infty}{lim}{x}_{i}(t)=0$ for

*i*∈[8,11]; - (iii)
If ${\mathcal{R}}_{v}>1$ and ${\mathcal{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 ${\mathcal{R}}_{v}$ and the Lyme disease invasion threshold ${\mathcal{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 ${\beta}_{{H}_{2}L}$ and ${\beta}_{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 ${\beta}_{{H}_{2}L}$ and ${\beta}_{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, ${\beta}_{{H}_{2}L}$=${\beta}_{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 ${\mu}_{{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.

*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, 47–49]

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*/*D*_{1}(*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 $\frac{1}{n}$ which is dependent of the temperatures of subsequent days.

*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]:

*θ*

^{ 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

*p*

_{1}=

*p*

_{2}=0.

## Results

We use various indices to measure the Lyme disease risk to humans: (i) ${\mathcal{R}}_{v}$, used to determine the tick population persistence; (ii) ${\mathcal{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 ${\mathcal{R}}_{0}<1$, while the seasonal risk pattern appears when ${\mathcal{R}}_{0}>1$. The numerical calculation of ${\mathcal{R}}_{v}$ is implemented by the dichotomy method where the system *d* *X*/*d* *t*=$(F(t)/{\mathcal{R}}_{v}-V(t))X(t)$ has a dominant Floquet multiplier equal to 1 [50]. A similar method is used to estimate ${\mathcal{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].

### 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.

*p*

_{2}=3.5) facilitate the growth of tick population and spread of the pathogen. For the Virginia opossum with a low reservoir competence (${\beta}_{{H}_{2}L}=0.0025$ and ${\beta}_{N{H}_{2}}=0.261$), the reduction of ${\mathcal{R}}_{0}$ largely attributes to not only the low transmission ability, but also their large biting biases coefficients (

*p*

_{1}=7.2,

*p*

_{2}=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 ${\mathcal{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.

*p*

_{1}and

*p*

_{2}. The result shows that ${\mathcal{R}}_{0}$ is very sensitive to the variations of both biting biases (Figure 7). Moreover, the relationships between ${\mathcal{R}}_{0}$ and

*p*

_{1}and

*p*

_{2}varies with host species: (i) ${\mathcal{R}}_{0}$ increases with increased

*p*

_{1}and

*p*

_{2}in the case of the eastern chipmunk, and therefore this species always facilitates disease transmission within our parameter region; (ii) the relation between ${\mathcal{R}}_{0}$ and the larvae bias

*p*

_{1}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.

## 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, ${\mathcal{R}}_{v}$ and ${\mathcal{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, 52–57] 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 (${\mathcal{R}}_{v}$, ${\mathcal{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 *x*^{0}∈*X*, system (1) admits a unique nonnegative solution *x*(*t*,*x*^{0}) through this initial value with the maximal interval of existence [0,*σ*) for some *σ*>0.

*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:

*f*(

*t*) with period

*τ*, denote $\hat{f}=\underset{t\in [0,\tau ]}{max}f(t)$ and $\stackrel{~}{f}=\underset{t\in [0,\tau ]}{min}f(t)$. It is easy to see that system (10) can be controlled by the following cooperative system:

If $\hat{{F}_{A}}\frac{\hat{{d}_{N}}}{{\mu}_{\mathit{\text{QA}}}}\frac{\hat{{F}_{N}}}{{\mu}_{\mathit{\text{FN}}}}\frac{\hat{{d}_{L}}}{{\mu}_{\mathit{\text{QN}}}}\frac{\hat{{F}_{L}}}{{\mu}_{\mathit{\text{FL}}}}\frac{\hat{{d}_{E}}}{{\mu}_{\mathit{\text{QL}}}}\frac{\hat{b}}{{\mu}_{E}}>{\mu}_{\mathit{\text{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*.

*τ*-positive periodic solution

*τ*-positive periodic solution (

*E*

^{∗}(

*t*), ${L}_{Q}^{\ast}(t)$, ${L}_{F}^{\ast}(t)$, ${N}_{Q}^{\ast}(t)$, ${N}_{F}^{\ast}(t)$, ${A}_{Q}^{\ast}(t)$, ${A}_{F}^{\ast}(t)$) is also

*τ*-periodic. Since for any $x\in {\mathbb{R}}_{+}^{7}\phantom{\rule{0.3em}{0ex}}\phantom{\rule{0.3em}{0ex}}\setminus \phantom{\rule{0.3em}{0ex}}\phantom{\rule{0.3em}{0ex}}\left\{0\right\}$, $\underset{\mathit{\text{n}}\to \infty}{lim}{P}^{6n}(x)$= (

*E*

^{∗}(0), ${L}_{Q}^{\ast}(0)$, ${L}_{F}^{\ast}(0)$, ${N}_{Q}^{\ast}(0)$, ${N}_{F}^{\ast}(0)$, ${A}_{Q}^{\ast}(0)$, ${A}_{F}^{\ast}(0))$ where

*P*is the Poincar$\stackrel{\u0301}{\mathrm{e}}$ map associated with the

*τ*-periodic system (3). Hence,

which implies that (*E*^{∗}(*t*), ${L}_{Q}^{\ast}(t)$, ${L}_{F}^{\ast}(t)$, ${N}_{Q}^{\ast}(t)$, ${N}_{F}^{\ast}(t)$, ${A}_{Q}^{\ast}(t)$, ${A}_{F}^{\ast}(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$\stackrel{\u0301}{\mathrm{e}}$ map of system (2), that is, *P*(*x*^{0})=*x*(11*τ*,*x*^{0}), where *x*(*t*,*x*^{0}) is the solution of system (2) through *x*^{0}. Then *P* is compact. Let *ω* = *ω*(*x*^{0}) be the omega limit set of *P*(*x*^{0}). 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*∈[1,9]. Hence,

*ω*= {(0, 0, 0, 0, 0, 0, 0, 0, 0)}×

*ω*

_{1}for some ${\omega}_{1}\subset {\mathbb{R}}^{2}$. It is easy to see that

*P*

_{1}is the Poincar$\stackrel{\u0301}{\mathrm{e}}$ map associated with the following equation:

Since *ω* is an internally chain transitive set for *P*, it easily follows that *ω*_{1} is an internally chain transitive set for *P*_{1}. 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.

*E*

^{∗}(

*t*), ${L}_{Q}^{\ast}(t)$, ${L}_{F}^{\ast}(t)$, ${N}_{Q}^{\ast}(t)$, ${N}_{F}^{\ast}(t)$, ${A}_{Q}^{\ast}(t)$, ${A}_{F}^{\ast}(t)$), for system (3) such that for any

*x*

^{0}with $\sum _{i=1}^{7}{x}_{i}^{0}>0$, we have

*P*

_{2}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

*P*

_{2}. Since ${\mathcal{R}}_{0}\le 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

Therefore, statement (ii) holds.

*x*

^{0}with $\sum _{i=1}^{7}{x}_{i}^{0}>0$, we have

where *P*_{2} 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 *P*_{2}. We claim that *ω*_{3}≠{0} for any $({x}_{8}^{0},{x}_{9}^{0},{x}_{10}^{0},{x}_{11}^{0})>0$.

*ω*

_{3}={0}. That is

*δ*>0 such that the spectral radius of the Poincar$\stackrel{\u0301}{\mathrm{e}}$ map associated with the following linearized system is greater than unity:

*u*

^{∗}(

*t*) such that

*τ*

_{0}>0 such that for all

*t*>

*τ*

_{0},

*t*>

*τ*

_{0}. By a standard comparison argument, we have

a contradiction to (13).

*ω*

_{3}≠{0} and the positive periodic solution $({L}_{\mathit{\text{FI}}}^{\ast}(t),{N}_{\mathit{\text{QI}}}^{\ast}(t),{H}_{1I}^{\ast}(t),{H}_{2I}^{\ast}(t))$ is globally asymptotically stable for system (5) in ${\mathbb{R}}_{+}^{4}\setminus \left\{0\right\}$, it follows that

*P*

_{2}. By [65], Theorem 3.1, we then get

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).