 Research
 Open Access
 Published:
Impact of biodiversity and seasonality on Lymepathogen transmission
Theoretical Biology and Medical Modelling volume 11, Article number: 50 (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 stagestructured 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 Lymepathogen transmission dynamics, are identified and shown to fully classify the longterm 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.
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 Lymepathogen 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 Lymepathogen 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 Lymepathogen 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 tickborne 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 tickborne 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) temperaturedependent/temperatureindependent development rates; (ii) temperaturedependent host seeking rates; (iii) densitydependent mortalities, caused by the hosts’ responses during the feeding period; (iv) constant mortalities of offhost 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.
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 stagestructured deterministic model is formulated in the next section. Moreover, two thresholds, one for the tick population dynamics and the other for the Lymepathogen 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 selfexplained as summarized in Table 1. For instance, L_{ FS } represents the subgroup of susceptible feeding larvae.
We assume that the host community of the tick population contains three species groups: (i) the whitefooted 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 whitetailed deer D (mainly Odocoileus virginianus), which is believed to be the paramount food provider for adults and intransmissible for the spread of Lymepathogen [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 Lymepathogen transmission cycle between the hosts and multistage tick population is presented in the diagram of Figure 1.
In the hostpathogentick 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, 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
Similarly, the nymphal infection rate which comes from the contact of questing susceptible nymphs and infectious hosts is given by
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
Similarly, the alternative host is infected by the infectious nymphal biting at a rate
Therefore, the disease transmission process between ticks and their hosts can be described by the following system:
We assume all the coefficients in the system are nonnegative and the timedependent coefficients are τperiodic with period τ=365 days. The detailed parameter definitions and sample values of these parameters are represented in Table 2. These timedependent parameters will be estimated in subsection ‘Timedependent parameters’ below.
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}\in X:=\{x\in {\mathbb{R}}_{+}^{14}:{H}_{1I}\le {H}_{1},{H}_{2I}\le {H}_{2}\}$, system (1) has a unique and bounded solution 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
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
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 stagestructured system for the tick population growth decoupled from system (2):
Linearization of system (3) at zero leads to the following linear system
Let 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) =
Then we can rewrite (4) as
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), 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
where 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.
Following ideas proposed in [39, 40], we define a next generation operator 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.
The following statements are valid:

(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
If threshold for ticks ${\mathcal{R}}_{v}>1$, then there exists a positive periodic solution,
for system (3) such that
In this case, equations for the infected populations in system (2) give rise to the following limiting system:
Following ideas of [39, 40], as proceed in the definition of ${\mathcal{R}}_{v}$ in the previous section, we can define a threshold for the pathogen. To do this, we introduce
and
Assume $\stackrel{~}{Y}$(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,
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 $\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 Lymepathogen 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.
Let 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 stagestructured 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 timedependent 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 shorttailed 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 noncompetent, ${\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 Lymepathogen since the species has a powerful immune system so that it can clean up the Lymepathogen 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.
Timedependent parameters
In order to investigate the impact of climate warming on the seasonal tick population abundance and Lymepathogen invasion, temperature is considered as a variable index in our study and we assume that it changes periodically with time. Therefore, those timedependent coefficients are indeed temperaturedependent, 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 timedependent 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 30year normal temperature data are collected from the Environment Canada website (Figure 2) [46].
Next we turn to the estimation of timedependent 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, 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 preoviposition 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.
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]:
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 timedependent parameters in oneyear 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 timedependent 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 timedependent 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 ${\mathcal{R}}_{v}$ increases from 1.38 to 1.62, and the values of ${\mathcal{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 ${\mathcal{R}}_{v}>1$ and establishment of pathogen population if ${\mathcal{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.
Impact of host biodiversity on the disease risk
Now, we seed the model with temperature condition in the 1981−2010 period so that timedependent 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 whitefooted mice and the whitetailed 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 densitydependent 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 ${\mathcal{R}}_{v}$ continuously increase with the increased number of hosts; while the change of ${\mathcal{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 ${\mathcal{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 (${\beta}_{{H}_{2}L}=0.569$ and ${\beta}_{N{H}_{2}}=0.971$), their ability of Lymepathogen transmission and high biting bias coefficient of nymphs (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 noncompetent host, but eventually reduces when the size of western fence lizard attains a certain level.
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 noncompetent 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 noncompetent host, does not always dilute the risk of Lyme disease.
We also perform sensitivity analysis of the threshold ${\mathcal{R}}_{0}$ against the biting biases 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 Lymepathogen, 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 reactiondiffusion 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 reactiondiffusion 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.
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:
For any periodic nonnegative function 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:
Clearly, there is only one nonnegative equilibrium zero for system (11) when
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.
Theorem 2.3.4 in [42] directly implies that if ${\mathcal{R}}_{v}\le 1$, then zero is globally asymptotically stable for system (3) in ${\mathbb{R}}_{+}^{7}$; if ${\mathcal{R}}_{v}>1$, then system (3) admits a unique 6τpositive periodic solution
and it is globally asymptotically stable for system (3) with initial values in ${\mathbb{R}}_{+}^{7}\phantom{\rule{0.3em}{0ex}}\phantom{\rule{0.3em}{0ex}}\setminus \left\{0\right\}$. It remains to prove that the 6τ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,
On the other hand,
Thus,
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) In the case where ${\mathcal{R}}_{v}\le 1$, we obtain $\underset{\mathit{\text{t}}\to \infty}{lim}{x}_{i}(t)=0$ for 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
where 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.
(ii) In the case where ${\mathcal{R}}_{v}>1$, then there exists a 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)$), for system (3) such that for any x ^{0} with $\sum _{i=1}^{7}{x}_{i}^{0}>0$, we have
Thus, $\omega =\{({E}^{\ast}(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))\}\times {\omega}_{2}$ for some ${\omega}_{2}\subset {\mathbb{R}}^{4}$, and
where 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.
(iii) In the case where ${\mathcal{R}}_{v}>1$ and ${\mathcal{R}}_{0}>1$, then there exists a 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))$, for system (3) such that for any x ^{0} with $\sum _{i=1}^{7}{x}_{i}^{0}>0$, we have
It then follows that $\omega =\{({E}^{\ast}(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)))\}\times {\omega}_{3}$ for some ${\omega}_{3}\subset {\mathbb{R}}^{4}$, and
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$.
Assume that, by contradiction, ω_{3}={0}. That is
for some $({x}_{8}^{0},{x}_{9}^{0},{x}_{10}^{0},{x}_{11}^{0})>0$. Then, we have
Since ${\mathcal{R}}_{0}>1$, there exists some δ>0 such that the spectral radius of the Poincar$\stackrel{\u0301}{\mathrm{e}}$ map associated with the following linearized system is greater than unity:
It then follows from the same argument as in the proof of Theorem 2.3 that the following system
admits a positive periodic u^{∗}(t) such that
Therefore, there exists some τ_{0}>0 such that for all t>τ_{0},
Hence, we conclude that