Mathematics in Applied Sciences and Engineering https://doi.org/10.5206/mase/10508 Volume 1, Number 1, March 2020, pp.65-84 https://ojs.lib.uwo.ca/mase MODELING THE IMPACT OF HOST RESISTANCE ON STRUCTURED TICK POPULATION DYNAMICS MAHNAZ ALAVINEJAD, JEMISA SADIKU, AND JIANHONG WU Abstract. For a variety of tick species, the resistance, behavioural and immunological response of hosts has been reported in the biological literature but its impact on tick population dynamics has not been mathematically formulated and analyzed using dynamical models reflecting the full biological stages of ticks. Here we develop and simulate a delay differential equation model, with a particular focus on resistance resulting in grooming behaviour. We calculate the basic reproduction number using the spectral analysis of delay differential equations with positive feedback, and establish the existence and uniqueness of a positive equilibrium when the basic reproduction number exceeds unit. We also conduct numerical and sensitivity analysis about the dependence of this positive equilibrium on the the parameter relevant to grooming behaviour. We numerically obtain the relationship between grooming behaviour and equilibrium value at different stages. 1. Introduction Lyme Disease is the most reported athropod-borne illness and it was first recognized in 1976 in Lyme, Connecticut USA [22]. Borrelia burgdorferi is a tick-borne spirochete responsible for Lyme disease which is found in nymphal Ixodes dammini and has the highest chance to be transmitted to the host if the infected tick feeds for a duration of 72 hours or more [16, 23, 10]. Once an infected tick bites the host, a skin lesion called erythema migrans (EM) starts emerging and more than 95% of those patients diagnosed with Lyme disease have EM on the tick biting site [5, 22]. Once the bacterium enters the body it starts spreading in many organs and tissues through the lymph system and blood [5]. As time progresses the patient will experience headache, neck pain, fever, fatigue, and migratory musculoskeletal pain [23, 10, 5]. The government of Canada has data representing an increase of 144 cases in 2009 to 2025 cases in 2017 [17]. The I.scopularis also known as a black-legged tick is the main carrier of B. burgdorferi and has a life cycle of nearly two years [22]. The tick population undergoes three main stages: L-larvae, N-nymph, A-adult and to move from one stage to the other ticks will quest feed and molt [13, 26, 18, 19]. Larvae and nymph feed on small rodents such as mice while adult ticks are more selective when it comes to their host since their body is larger compared to larvae and nymph and therefore the host must be a large mammal such as a deer. For ticks to move from one stage to the next it requires three hosts per stage and often the tick may use the same host for all three blood meals [13, 26, 18, 19]. Female ticks lay eggs in the spring and larvae hatch during late summer. The larvae that feeds during the late summer starts molting to nymph during winter. The nymph then Received by the editors 12 March 2020; revised 29 March 2020; accepted 29 March 2020; published online 31 March 2020. 2000 Mathematics Subject Classification. Primary 92D25; Secondary 34K20. Key words and phrases. Ticks, host resistance, population dynamics. This research was supported by Natural Sciences and Engineering Research Council of Canada and the Canada Research Chair Program. 65 66 M. ALAVINEJAD, J. SADIKU, AND J. WU starts feeding in the spring of the following year and molts into adult on the same year. Adult ticks die shortly right after they lay their eggs in the early spring [26, 18]. When a tick bites a host the expression of immunity varies depending on different hosts and tick species.The effects on ticks can vary from a simple rejection of the tick to interfering with the duration of feeding, inhibition of egg laying, also decreasing their viability to death of the tick while feeding. In addition, studies reveal that when female ticks feed on immune cattle their body of fully engorged tick was reduced by 30% [12, 24, 2]. According to Brown [3] hosts with resistance respond to tick bites with an intensified grooming behaviour and the attachment site is marked by serous exudes which could engulf the tick. In an experiment conducted on resistant guinea pigs bitten by Dermacentor andersoni, basophilia is present on the biting site. The attachment of a tick on a tick-sensitizes host is characterized by packs of basophils located in the intraepidermal vesicles. When ticks’ extracts are injected into tick- sensitized host it causes a skin reaction and the plasma of the host expresses anti-tick antibodies which suggests a present mediated immune response. In case of unbitten animals, the reaction starts with neutrophils and the feeding site is characterized by an hemorrhagic as feeding progresses. Basophils start to also accumulate to the feeding site, however little degranulation occurs. In an experiment to study the effect of resistance of guinea pigs to ticks, basophil degranulation at tick feeding sites, resulted in tick rejection after tick-attachment: 29% after 6 hours, 18% after 12 hours, 22% after 24 hours, 37% after 48 hours and 7.3% after 72 and 96 hours. This shows that ticks are most susceptible to the resistance at 6, 12, 24 and 48 hours after attachment which are corresponding to the attachment time and fast feeding period [4]. There have been intensive studies modelling the dynamics of tick-host interaction and the transmis- sion of various pathogens. Different aspects have also been included such as: seasonality , environmental changes, geographical heterogeneity and so on. On the other hand, few models incorporate delays in the development of tick from each life stage to the next [6, 25, 27]. Jennings et al. [9] studied the effect of host resistance on tick population dynamics. They developed a mathematical model, described by a system of ordinary differential equations, focusing on tick-host interaction where the tick’s life cycle was divided into two main stages, adult and juvenile, and the host was subdivided into host with no immunity and host with immunity. Their focus is to show how immunity affects the extinction or per- sistence of tick dynamics. However, their model does not include all biological stages (and sub-stages) of ticks and the possibility of different immunological response for each stage is ignored. Here, we consider a stage-structured model that involves the full biological dynamics of tick and the emphasis is on the grooming behaviour of the host and its impact on tick population dynamics. We analyze the grooming behaviour in the mathematical model as a reduction in the successful attachment rates of ticks on the host i.e., the host-finding rates are reduced by a fraction for the host that shows resistance to tick bites. The model studies three main stages of tick’s life cycle in which the ticks interact with hosts during questing-feeding-molting process. There is one more stage that we consider between Adult and Egg which is egg laying female. The host is divided into two compartments: host with resistance (host has been bitten by ticks before) and host with no resistance (host that has not been exposed to ticks). We observe that the basic reproduction number does not change with the resistance factor, however, numerical simulations show that the value of the positive equilibrium for different stages of tick population, and the dynamical behaviour of the solutions change with varying the resistance factor. Also, the sensitivity analysis demonstrates the dependence of the solutions on different parameters. MODELING THE IMPACT OF HOST RESISTANCE ON STRUCTURED TICK POPULATION DYNAMICS 67 Lq Lf Lm Nq Nf Nm AqAfAelfE βL(αLHr+ + Hr−)Lq δLf δe−d Lmτ(L,N)Lf (t− τ(L,N)) βN (αNHr+ + Hr−)Nq δLf δe−d Nmτ(N,A)Nf (t− τ(N,A)) βA(αAHr+ + Hr−)AqεδAf γ(Aelf ) γ (A e lf (t − τ ( E ,L ) )) d L q L q d L m L m d N q N q d N m N m d A q A q d A e l f A e lfd E E Figure 1. Flow diagram for the ticks’ life cycle and their interaction with hosts 2. The Model Formulation We start the model aiming to focus on the grooming behaviour. We model the three stages of larvae, nymph, and adult. The larvae and nymph populations are subdivided into questing, feeding and molting. On the other hand, for the adult population we consider adult egg laying female Aelf , adult questing (Aq) and feeding (Af ). Since a single female tick lays several thousands eggs the birth rate is the entry into population which is represented by Ricker function, γ(A) = pAe−qA [19, 15]. Tick dynamics are regulated by insufficient resources for blood meal and this is illustrated in parameter q. The delay functions, demonstrating the time delays of ticks molting from one stage to another, are constants. In the model, τ(E,L), τ(L,N), τ(N,A) represent the time it takes for ticks to molt from egg to larvae, larvae to nymph and nymph to adult, respectively. The host population is divided into two compartments: Hr+ represents the bitten host compartment that have developed with immunity; Hr− represents the compartment of hosts that have never been bitten before and therefore without immunity. H is the total host population with a birth rate b and a mortality rate µ. The density- dependent regulations of the host population is described by K, c, and b − µ. The variables and parameters and their meaning are given in Tables 3 and 1. The life cycle of ticks and their interaction with hosts is illustrated in Figure 1. We suppose the successful attachment rates are reduced by a fraction αL for larvea, αN for nymph and αA for adult ticks. Based on the results from [4] we can assume that α is in the range [0.6, 0.95], however we will study the effect of α values in [0, 1]. We also assume the hosts with no resistance develop resistance to ticks at a rate, denoted by κ, that depends on the tick densities, tick attachment rates and the immune system response. 68 M. ALAVINEJAD, J. SADIKU, AND J. WU Table 1. Definition of parameters and their values Symbol Meaning Value Reference dLq Per capita mortality rate of Lq 0.6 × 10−2 per day [15] dLm Per capita mortality rate of Lm 0.3 × 10−2 per day [15] dNq Per capita mortality rate of Nq 0.6 × 10−2 per day [15] dNm Per capita mortality rate of Nm 0.2 × 10−2 per day [15] dAq Per capita mortality rate of Aq 0.6 × 10−2 per day [15] dAelf Per capita mortality rate of Aelf 1 per day [28] dE Per capita mortality rate of E 0.2 × 10−2 per day [15] βL Successful attachment rate of 0.6 × 10−3 per day per host [11] questing larva to host βN Successful attachment rate of 0.6 × 10−3 per day per host [11] questing nymph to host βA Successful attachment rate of 0.2 × 10−2 per day per host [11] questing adult to host β∗L Rate of developing resistance to larva κ×βL per day per tick Calculated β∗N Rate of developing resistance to nymph κ×βN per day per tick Calculated β∗A Rate of developing resistance to adult κ×βA per day per tick Calculated δ Detachment rate 0.01 per day [20] αL Host grooming effect for larva 0.4, [0, 1] unitless Assumed αN Host grooming effect for nymph 0.6, [0, 1] unitless Assumed αA Host grooming effect for adult 0.5, [0, 1] unitless Assumed � Female proportion 0.5 unitless [7] τ(E,L) The delay of development 21 days [15] form egg to larvae τ(L,N) The delay of development 101.18 ×Temp−2.25, 200 days [15] form larvae to nymph τ(N,A) The delay of development 1596 ×Temp−1.21, 61 days [15] form nymph to adult b Birth rate of the host 0.66 × 10−3 per day [25] µ Death rate of the host 0.33 × 10−3 per day [25] c Crowding 3.5 × 10−4 per day Calculated K Carrying Capacity of deers 20 [15] p Maximum number of eggs 3000 [15] per female adult tick q The strength of density dependence 0.001 unitless [19] κ Constant factor for resistance development 0.0001 unitless Assumed MODELING THE IMPACT OF HOST RESISTANCE ON STRUCTURED TICK POPULATION DYNAMICS 69 Table 2. Modified parameter values to get different values for R0 Symbol Modified Value Comments dE 1.2 × 0.2 × 10−2 p+20%p dLq 1.2 × 0.6 × 10−2 p+20%p dLm 1.2 × 0.3 × 10−2 p+20%p dNq 1.2 × 0.6 × 10−2 p+20%p dNm 1.2 × 0.2 × 10−2 p+20%p dAq 1.2 × 0.6 × 10−2 p+20%p βL 0.1 × 0.6 × 10−3, 0.2 × 0.6 × 10−3 10%p, 20%p βN 0.3 × 0.6 × 10−3, 0.5 × 0.6 × 10−3 30%p, 50%p βA 0.5 × 0.2 × 10−2 50%p fixed β∗L κ×βL changed by changing βL β∗N κ×βN changed by changing βN β∗A κ×βA changed by changing βA αL 0.4 varied in [0, 1] αN 0.6 varied in [0, 1] αA 0.5 varied in [0, 1] c 1.2 × 3.5 × 10−4 p + 20%p fixed q 0.001 not changed Table 3. Definition of Variables and their initial values Symbol Meaning Initial value Lq Number of questing larvae Lq0 = 1 × 106 Lf Number of feeding larvae Lf0 (θ) = 0, −τ(E,L) ≤ θ ≤ 0 Lm Number of molting larvae Nq Number of questing nymph Nq0 = 0 Nf Number of feeding nymph Nf0 (θ) = 0, −τ(L,N) ≤ θ ≤ 0 Nm Number of molting nymph Aq Number of questing adult Aq0 = 0 Af Number of feeding adult Af0 = 0 Aelf Number of egg laying female adult Aelf0 (θ) = 0, −τ(N,A) ≤ θ ≤ 0 E Number of eggs H Number of hosts Hr+ Number of hosts with resistance Hr+ = 0 Hr− Number of host with no resistance 70 M. ALAVINEJAD, J. SADIKU, AND J. WU In order to make the model comprehensible we neglect few biological factors of tick dynamics. There are multiple blood meals that take place during molting procedures however in our model we consider only a homogeneous molting process, that is, ticks feed once, drop and molt with a constant time delay. The death rate depends on the stage of the tick (egg, larvae, nymph, adult) and also on whether the tick is questing or feeding. However, we consider a constant mortality rate. Impact of climate change on development of ticks having a non linear relationship with increasing ambient temperature has not also been modelled. In addition, the questing rate is considered constant, even though it decreases as the temperatures and the day light decreases. The model is described by the following system of delay differential equations:   dLq dt = e−d Eτ(E,L)γ(Aelf (t− τ(E,L))) −βLLq(t)(αLHr+(t) + Hr−(t)) −dLqLq(t) dLf dt = βLLq(t)(αLHr+(t) + Hr−(t)) − δLf (t) dLm dt = δLf (t) −dLmLm(t) − δψe−d Lmτ(L,N)Lf (t− τ(L,N)) dNq dt = δψe−d Lmτ(L,N)Lf (t− τ(L,N)) −βNNq(t)(αNHr+(t) + Hr−(t)) −dNqNq(t) dNf dt = βNNq(t)(αNHr+(t) + Hr−(t)) − δNf (t) dNm dt = δNf (t) −dNmNm(t) − δψe−d Nmτ(N,A)Nf (t− τ(N,A)) dAq dt = δψe−d Nmτ(N,A)Nf (t− τ(N,A)) −βAAq(t)(αAHr+(t) + Hr−(t)) −dAqAq(t) dAf dt = βAAq(t)(αAHr+(t) + Hr−(t)) −δAf (t) dAelf dt = εδAf (t) −dAelf Aelf (t) dE dt = γ(Aelf (t)) −dEE(t) −e−d Eτ(E,L)γ(Aelf (t− τ(E,L))) dHr− dt = bH(t) −µHr−(t) − c K H(t)Hr−(t) − (β∗LLq(t) + β ∗ NNq(t) + β ∗ AAq(t))Hr−(t) dHr+ dt = −µHr+(t) − c K H(t)Hr+(t) + (β ∗ LLq(t) + β ∗ NNq(t) + β ∗ AAq(t))Hr−(t) (2.1) where γ(A) = pAe−qA is the birth function. Here, we use the following equation for the host population dynamics dH(t) dt = (b−µ)H(t) − c K (H(t))2 (2.2) where H(t) = Hr−(t) + Hr+(t). Note that the positive equilibrium of this equation is given by H̄ = (b−µ) c K. Interpreting K as an environmental constraint, and in order to have H̄ ≤ K we assume c ≥ (b−µ), with H̄ = K when the equality holds. MODELING THE IMPACT OF HOST RESISTANCE ON STRUCTURED TICK POPULATION DYNAMICS 71 From System (2.1) we can get the following integral equations for Lm(t), Nm(t) and E(t) Lm(t) = Lm(0) − ∫ 0 −τ(L,N) e−d Lm (−s)δLf (s)ds + ∫ t t−τ(L,N) e−d Lm (t−s)δLf (s)ds Nm(t) = Nm(0) − ∫ 0 −τ(N,A) e−d Nm (−s)δNf (s)ds + ∫ t t−τ(N,A) e−d Nm (t−s)δNf (s)ds E(t) = E(0) − ∫ 0 −τ(E,L) e−d E (−s)γ(Aelf (s))ds + ∫ t t−τ(E,L) e−d E (t−s)γ(Aelf (s))ds (2.3) Therefore System (2.1) is equivalent to the following  dLq dt = e−d Eτ(E,L)γ(Aelf (t− τ(E,L))) −βLLq(t)(αLHr+(t) + Hr−(t)) −dLqLq(t) dLf dt = βLLq(t)(αLHr+(t) + Hr−(t)) − δLf (t) dNq dt = δψe−d Lmτ(L,N)Lf (t− τ(L,N)) −βNNq(t)(αNHr+(t) + Hr−(t)) −dNqNq(t) dNf dt = βNNq(t)(αNHr+(t) + Hr−(t)) − δNf (t) dAq dt = δψe−d Nmτ(N,A)Nf (t− τ(N,A)) −βAAq(t)(αAHr+(t) + Hr−(t)) −dAqAq(t) dAf dt = βAAq(t)(αAHr+(t) + Hr−(t)) −δAf (t) dAelf dt = εδAf (t) −dAelf Aelf (t) dHr− dt = bH(t) −µHr−(t) − c K H(t)Hr−(t) − (β∗LLq(t) + β ∗ NNq(t) + β ∗ AAq(t))Hr−(t) dHr+ dt = −µHr+(t) − c K H(t)Hr+(t) + (β ∗ LLq(t) + β ∗ NNq(t) + β ∗ AAq(t))Hr−(t) (2.4) together with (2.3). For further analyses of this model we use the theory of monotone dynamical systems [21]. Let τ = (τ1, · · · ,τ12) where τi ≥ 0, τ2 = τ(L,N), τ5 = τ(N,A)), τ9 = τ(E,L) are non zero and τi = 0 for i 6= 2, 5, 9. Assume |τ| = max{τi}. Let Cτ be the product of Banach spaces Cτi = C([−τi, 0],R), i.e., Cτ = 12∏ i=1 C([−τi, 0],R). Let Xt = (X 1 t , · · · ,X12t ) ∈ Cτ be given by Xit(θ) = X i(t + θ), i = 1, · · · , 12. where X(t) = (X1(t), · · · ,X12(t)) = (Lq,Lf,Lm,Nq,Nf,Nm,Aq,Af,Aelf,E,Hr−,Hr+). Then the right hand side of the Equation (2.1) is given by X′(t) = f(Xt). (2.5) We assume the initial data is non-negative. So we will assume the initial data X0 is in the Banach space C+τ defined below C+τ = {φ ∈ Cτ : φi(θ) ≥ 0,−τi ≤ θ ≤ 0}. 72 M. ALAVINEJAD, J. SADIKU, AND J. WU Also, for the initial data to be continuous and positive we assume: Lm(0) ≥ ∫ 0 −τ(L,N) e−d Lm (−s)δLf (s)ds Nm(0) ≥ ∫ 0 −τ(N,A) e−d Nm (−s)δNf (s)ds E(0) ≥ ∫ 0 −τ(E,L) e−d E (−s)γ(Aelf (s))ds. (2.6) The fundamental theory of functional differential equations implies that the solutions exist and are unique for all t ≥ 0. We now show that the solutions will be positive and remain bounded. Theorem 2.1. Let Xi(0) > 0 and Xi(θ) ≥ 0 for −τi ≤ θ < 0, for i = 1, · · · , 12. Then the solutions to the System (2.4) are positive and bounded for all t ≥ 0. Proof. Consider the first equation in (2.4). First we look at the solution on [0,τ]: if there exists t1 ∈ (0,τ) such that Lq(t1) = 0, then the derivative dLq(t)/dt at t1 is dLq(t) dt ∣∣∣∣ t1 = e−d Eτ(E,L)γ(Aelf (t1 − τ(E,L))). (2.7) Since initial data for Aelf on [−τ, 0] is positive, the derivative of Lq at t1 is positive and therefore Lq(t) is increasing, so it can not be negative. The same argument can be applied for [τ, 2τ]. This proves that Lq(t) ≥ 0 for all t ≥ 0. If there exists t2 such that Lf (t2) = 0, then the derivative of Lf at t2 is given by dLf (t) dt ∣∣∣∣ t2 = βLLq(t)(αLHr+(t) + Hr−(t)) (2.8) which is positive since Lq(t2), Hr+(t2) and Hr−(t) are positive. Thus Lf is increasing at t2 so it can not be negative. The same argument applies for other equations. Therefore the solutions are positive. From Equation (2.2) it is clear that H(t) is positive and bounded by the carrying capacity K. Also the above discussion shows that Hr− and Hr+ are positive for all t ≥ 0. We show the boundedness of the tick population as follows. Let T > 0 and τ = max{τ(E,L),τ(L,N),τ(N,A)}. We integrate the first equation in the original system (2.1) Lq(t) = e −dLq t−βL ∫ t 0 (αLHr+(u)+Hr−(u))du ∫ t 0 ed Lq s+βL ∫ s 0 (αLHr+(u)+Hr−(u))du ( e−d Eτ(E,L)γ(Aelf (s− τ(E,L))) ) ds + e−d Lq t−βL ∫ t 0 (αLHr+(u)+Hr−(u))duLq(0) therefore sup −τ≤t≤T Lq(t) ≤ Lq(0) + e−d Eτ(E,L) dLq sup −τ≤t≤T γ(Aelf (t)) using sup −τ≤t≤T e−βL ∫ t s (αLHr+(u)+Hr−(u))du = 1 and ∫ t 0 e−d Lq (t−s)ds < 1/dLq. Using the fact that γ(Aelf (t)) ≤ p/qe for all t ≥ 0, we see that sup −τ≤t≤T Lq(t) ≤ C where C = Lq(0) + pe −dEτ(E,L)/qedLq is independent of T . Therefore Lq(t) ≤ C for all −τ ≤ t < ∞. MODELING THE IMPACT OF HOST RESISTANCE ON STRUCTURED TICK POPULATION DYNAMICS 73 Integrating the next equations and taking the supermom we have: sup −τ≤t≤T Lf (t) ≤ sup −τ≤t≤0 Lf0 (t) + βLK δ sup −τ≤t≤T Lq(t) sup −τ≤t≤T Nq(t) ≤ Nq(0) + e−d Lmτ(L,N) dNq sup −τ≤t≤T Lf (t) sup −τ≤t≤T Nf (t) ≤ sup −τ≤t≤0 Nf0 (t) + βNK δ sup −τ≤t≤T Nq(t) sup −τ≤t≤T Aq(t) ≤ Aq(0) + e−d Nmτ(N,A) dAq sup −τ≤t≤T Nf (t) sup −τ≤t≤T Af (t) ≤ Af (0) + βAK δ sup −τ≤t≤T Aq(t) sup −τ≤t≤T Aelf (t) ≤ sup −τ≤t≤0 Aelf0 (t) + �δ dAelf sup −τ≤t≤T Af (t). Combining the above inequalities and assuming that the initial data are bounded we can see that these tick stages are bounded on −τ ≤ t < ∞. We can get similar inequalities from System (2.3). This proves that all tick stages are bounded. � Since the host population stabilizes quickly at H̄ = (b−µ)K/c, the limiting system is as follows   dLq dt = e−d Eτ(E,L)γ(Aelf (t− τ(E,L))) + βL(1 −αL)Lq(t)Hr+(t) − (βLH̄ + dLq )Lq(t) dLf dt = −βL(1 −αL)Lq(t)Hr+(t) + βLH̄Lq(t) −δLf (t) dNq dt = δψe−d Lmτ(L,N)Lf (t− τ(L,N)) + βN (1 −αN )Nq(t)Hr+(t) − (βNH̄ + dNq )Nq(t) dNf dt = −βN (1 −αN )Nq(t)Hr+(t) + βNH̄Nq(t) −δNf (t) dAq dt = δψe−d Nmτ(N,A)Nf (t− τ(N,A)) + βA(1 −αA)Aq(t)Hr+(t) − (βAH̄ + dAq )Aq(t) dAf dt = −βAAq(t)(1 −αA)Hr+(t) + βAH̄Aq(t) − δAf (t) dAelf dt = εδAf (t) −dAelf Aelf (t) dHr+ dt = −µHr+(t) − c K H̄Hr+(t) + (β ∗ LLq(t) + β ∗ NNq(t) + β ∗ AAq(t))(H̄ −Hr+(t)) (2.9) From now on we refer to this system as the main system of our model unless otherwise stated. 3. Analyses In this section we give the necessary condition for existence and uniqueness of the positive equilibrium point and the conditions for local stability of the tick free equilibrium. 74 M. ALAVINEJAD, J. SADIKU, AND J. WU 3.1. Equilibria. Let X∗ denote the vector (Lq,Lf,Nq,Nf,Aq,Af,Aelf,Hr+) in R8, and let f(X) be the right hand side of (2.9). In order to find all equilibria we need to solve the system f(X) = 0:  0 = e−d Eτ(E,L)γ(Aelf (t− τ(E,L))) + βL(1 −αL)Lq(t)Hr+(t) − (βLH̄ + dLq )Lq(t) 0 = −βL(1 −αL)Lq(t)Hr+(t) + βLH̄Lq(t) − δLf (t) 0 = δψe−d Lmτ(L,N)Lf (t− τ(L,N)) + βN (1 −αN )Nq(t)Hr+(t) − (βNH̄ + dNq )Nq(t) 0 = −βN (1 −αN )Nq(t)Hr+(t) + βNH̄Nq(t) − δNf (t) 0 = δψe−d Nmτ(N,A)Nf (t− τ(N,A)) + βA(1 −αA)Aq(t)Hr+(t) − (βAH̄ + dAq )Aq(t) 0 = −βAAq(t)(1 −αA)Hr+(t) + βAH̄Aq(t) − δAf (t) 0 = εδAf (t) −dAelf Aelf (t) 0 = −µHr+(t) − c K H̄Hr+(t) + (β ∗ LLq(t) + β ∗ NNq(t) + β ∗ AAq(t))(H̄ −Hr+(t)) (3.1) At the tick-free equilibrium, where all tick stages are equal to zero, we have Hr+ = 0. Let Hr+ 6= H̄ so that (H̄−(1−αL)Hr+), (H̄−(1−αN )Hr+), (H̄−(1−αA)Hr+) > 0. We want to derive conditions for existence and uniqueness of a (strongly) positive equilibrium point (Xi > 0 for all i = 1, · · · ,n). From the equations in (3.1) we get the following Lq = dAelf (βN (H̄ − (1 −αN )Hr+) + dNq )(βA(H̄ − (1 −αA)Hr+) + dAq ) s2s3�βLβNβL(H̄ − (1 −αL)Hr+)(H̄ − (1 −αN )Hr+)(H̄ − (1 −αA)Hr+) Aelf Lf = dAelf (βN (H̄ − (1 −αN )Hr+) + dNq )(βA(H̄ − (1 −αA)Hr+) + dAq ) s2s3δ�βNβA(H̄ − (1 −αN )Hr+)(H̄ − (1 −αA)Hr+) Aelf Nq = dAelf (βA(H̄ − (1 −αA)Hr+) + dAq ) s3�βNβA(H̄ − (1 −αN )Hr+)(H̄ − (1 −αA)Hr+) Aelf Nf = dAelf (βA(H̄ − (1 −αA)Hr+) + dAq ) s3δ�βA(H̄ − (1 −αA)Hr+) Aelf Aq = dAelf �βA(H̄ − (1 −αA)Hr+) Aelf Af = dAelf �δ Aelf (3.2) where s1 = e −dEτ(E,L) , s2 = ψe −dLmτ(L,N) and s3 = ψe −dNmτ(N,A) . From the first equation in the system (3.1) we get Lq = s1γ(Aelf ) (βL(H̄ − (1 −αL)Hr+) + dLq ) (3.3) and therefore γ(Aelf ) = d Aelf (βL(H̄−(1−αL)Hr+)+dLq )(βN (H̄−(1−αN )Hr+)+dNq )(βA(H̄−(1−αA)Hr+)+dAq ) s1s2s3�βLβNβA(H̄−(1−αL)Hr+)(H̄−(1−αN )Hr+)(H̄−(1−αA)Hr+) Aelf . (3.4) MODELING THE IMPACT OF HOST RESISTANCE ON STRUCTURED TICK POPULATION DYNAMICS 75 Since γ(Aelf ) = pAelfe −qAelf we have the following cases: Aelf = 0 or pe−qAelf = d Aelf (βL(H̄−(1−αL)Hr+)+dLq )(βN (H̄−(1−αN )Hr+)+dNq )(βA(H̄−(1−αA)Hr+)+dAq ) s1s2s3�βLβNβA(H̄−(1−αL)Hr+)(H̄−(1−αN )Hr+)(H̄−(1−αA)Hr+) (3.5) Finally, we reduce the system (3.1) to the following system 0 = Γ(Hr+) −pe−qAelf (3.6a) 0 = −bHr+ + Ω(Hr+)(H̄ −Hr+)Aelf (3.6b) where Γ(Hr+) = dAelf (βL(H̄ − (1 −αL)Hr+) + dLq )(βN (H̄ − (1 −αN )Hr+) + dNq )(βA(H̄ − (1 −αA)Hr+) + dAq ) s1s2s3�βLβNβA(H̄ − (1 −αL)Hr+)(H̄ − (1 −αN )Hr+)(H̄ − (1 −αA)Hr+) Ω(Hr+) = β ∗ L dAelf (βN (H̄ − (1 −αN )Hr+) + dNq )(βA(H̄ − (1 −αA)Hr+) + dAq ) s2s3�βLβNβA(H̄ − (1 −αL)Hr+)(H̄ − (1 −αN )Hr+)(H̄ − (1 −αA)Hr+) + β∗N dAelf (βA(H̄ − (1 −αL)Hr+) + dAq ) s3�βNβA(H̄ − (1 −αN )Hr+)(H̄ − (1 −αA)Hr+) + β∗A dAelf �βA(H̄ − (1 −αA)Hr+) From (3.6b) we have Aelf = bHr+ Ω(Hr+)(H̄ −Hr+) given that Hr+ 6= H̄ and Ω(Hr+) 6= 0 (it can be proved that this holds). Substituting this in the equation (3.6a) we get the following Γ(Hr+) = pe −q bHr+ Ω(Hr+)(H̄−Hr+). (3.7) This is a nonlinear equation for Hr+and we need to determine under what conditions this equation has a unique positive solution. Let G(Hr+) be the right hand side of Equation (3.7). The functions Γ and G have the following properties: (i) Γ is a rational function and is strictly increasing for 0 < Hr+ < H̄; (ii) Γ(0) > 0 is given by dAelf (βLH̄ + d Lq )(βAH̄ + d Aq )(βNH̄ + d Nq ) s1s2s3�βLβNβAH̄3 ; (iii) G is a negative exponential function and it approaches zero (exponentially) as Hr+ approaches H̄; (iv) G(0) = p. From these properties we can see that the equation (3.7) has at least one solution 0 < Hr+ < H̄, if and only if G(0) > Γ(0), i.e., p > dAelf (βLH̄ + d Lq )(βAH̄ + d Aq )(βNH̄ + d Nq ) s1s2s3�βLβNβAH̄3 . This solution is unique if G(Hr+) is monotonically decreasing, and this holds if and only if d dHr+ ( Hr+ Ω(Hr+)(H̄ −Hr+) ) > 0 for all Hr+ ∈ (0,H̄). 76 M. ALAVINEJAD, J. SADIKU, AND J. WU Theorem 3.1. Let Rv0 = ps1s2s3�βLβNβAH̄ 3 dAelf (βLH̄ + d Lq )(βAH̄ + d Aq )(βNH̄ + d Nq ) . If Rv0 > 1, then system (2.9) has a positive equilibrium point. If additionally d dHr+ ( Hr+ Ω(Hr+)(H̄ −Hr+) ) > 0 holds, then the positive equilibrium is unique. 3.2. Stability of the tick-free Equilibrium. First we linearize System (2.9) about a given equilib- rium point using the Fréchet derivative of the function F(X), given by the right hand side of the System (2.9): DF(X∗)X = lim h→0 (F(X∗ + hX) −F(X∗) h ) The linearized system is given by Df1(X ∗)X = pe−d Eτ(E,L)Aelf (t− τ(E,L))e−qA ∗ elf (t−τ(E,L)) −pqe−d Eτ(E,L)Aelf (t− τ(E,L))A∗elf (t− τ(E,L))e −qA∗elf (t−τ(E,L)) + (1 −αL)βL(L∗q(t)Hr+(t) + Lq(t)H ∗ r+(t)) − (βLH̄ + d Lq )Lq(t) Df2(X ∗)X = −(1 −αL)βL(L∗q(t)Hr+(t) + Lq(t)H ∗ r+(t)) + βLH̄Lq(t) −δLf (t) Df4(X ∗)X = δψe−d Lmτ(L,N)Lf (t− τ(L,N)) + (1 −αN )βN (N∗q (t)Hr+(t) + Nq(t)H ∗ r+(t)) − (βNH̄ + d Nq )Nq(t) Df5(X ∗)X = −(1 −αN )βN (N∗q (t)Hr+(t) + Nq(t)H ∗ r+(t)) + βNH̄Nq(t) − δNf (t) Df7(X ∗)X = δψe−d Nmτ(N,A)Nf (t− τ(N,A)) + (1 −αA)βA(A∗q(t)Hr+(t) + Aq(t)H ∗ r+(t)) − (βAH̄ + d Aq )Aq(t) Df8(X ∗)X = −(1 −αA)βA(A∗q(t)Hr+(t) + Aq(t)H ∗ r+(t)) + βAH̄Aq(t) −δAf (t) Df9(X ∗)X = εδAf (t) −dAelf Aelf (t) Df12(X ∗)X = −(µ + c K H̄)Hr+(t) + H̄(β ∗ LL ∗ q(t) + β ∗ NN ∗ q (t) + β ∗ AA ∗ q(t)) − ( (β∗LLq(t) + β ∗ NNq(t) + β ∗ AAq(t))H ∗ r+(t) + (β ∗ LL ∗ q(t) + β ∗ NN ∗ q (t) + β ∗ AA ∗ q(t))Hr+(t) ) (3.8) The linearized system about the tick-free equilibrium point is as follows:  Df1(X ∗)X = ps1Aelf (t− τ(E,L)) − (βLH̄ + dLq )Lq(t) Df2(X ∗)X = βLH̄Lq(t) − δLf (t) Df4(X ∗)X = δs2Lf (t− τ(L,N)) − (βNH̄ + dNq )Nq(t) Df5(X ∗)X = βNH̄Nq(t) − δNf (t) Df7(X ∗)X = δs3Nf (t− τ(N,A)) − (βAH̄ + dAq )Aq(t) Df8(X ∗)X = βAH̄Aq(t) −δAf (t) Df9(X ∗)X = εδAf (t) −dAelf Aelf (t) Df12(X ∗)X = −(µ + c K H̄)Hr+(t) (3.9) Using the theory of monotone dynamical systems we can see that system (2.9) is cooperative ([21] corollary 3.2) and therefore stability of the zero equilibrium of system (3.9) is given by the stability of the corresponding ODE system. MODELING THE IMPACT OF HOST RESISTANCE ON STRUCTURED TICK POPULATION DYNAMICS 77 Theorem 3.2. If Rv0 < 1, then X = 0 is the only equilibrium point of the system (2.9) and is locally asymptotically stable. When Rv0 > 1, there exists a positive equilibrium point and X = 0 is unstable. Proof. We use the method of next generation matrix for the ODE system given by X′(t) = JX(t) where the matrix J is obtained from system (3.9): J =   −βLH̄ −dLq 0 0 0 0 0 ps1 0 βLH̄ −δ 0 0 0 0 0 0 0 δs2 −βNH̄ −dNq 0 0 0 0 0 0 0 βNH̄ −δ 0 0 0 0 0 0 0 δs3 −βAH̄ −dAq 0 0 0 0 0 0 0 βAH̄ −δ 0 0 0 0 0 0 0 �δ −dAelf 0 0 0 0 0 0 0 0 −b   The matrix J can be written as J = F − V . The zero equilibrium is locally asymptotically stable if ρ(FV −1) < 1 (ρ is the spectral radius of FV −1) and it is unstable if ρ(FV −1) > 1. We can see that ρ(FV −1) = Rv0 = ps1s2s3�βLβNβAH̄ 3 dAelf (βLH̄ + d Lq )(βAH̄ + d Aq )(βNH̄ + d Nq ) . � 4. Numerical Simulations In this section we study the long-term dynamical behaviour of the system using numerical simulations and perform a sensitivity analysis for different parameters. 4.1. Model parametrization and validation. The observation of the dynamical behaviour of each stage of the tick population is demonstrated by applying DDE23 packages in Matlab to System (2.9). The model is parameterized using parameter values available in mathematical and ecological literature ([7, 11, 15, 20, 19, 28]). Parameter values and initial conditions are given in Tables 1-3. We note that the grooming behaviour does not impact the initial growth of the tick population, since parameters reflecting the grooming factor do not change the value of the basic reproduction number. We consider three cases to illustrate the dynamics of tick population in the presence of grooming factor. However, in these cases we fix the values for parameters related to the grooming behaviour. In the first case (Figure 2) the basic reproduction number is below the threshold value i.e., Rv0 < 1, the tick-free equilibrium is locally asymptotically stable and therefore all stages of ticks go extinct. In case 2 (Figure 3) the basic reproduction number is slightly greater than one, the tick-free equilibrium point becomes unstable and the solutions approach the positive equilibrium without any initial oscillatory behaviour. In case 3 (Figure 5) the solutions oscillate initially and then approach the positive equilibrium. When the resistance related parameter values are fixed and the rest of the parameters vary, the positive equilibrium becomes unstable and a limit cycle appears. Therefore, the solutions oscillate about the equilibrium point. Figure 4 shows how a limit cycle appears as the value of αA increases from 0 to 1. To study the population behaviour without grooming factor we set αL = αN = αA = 1 and κ = 0 and for intense grooming behaviour the αL = αN = αA = 0. In addition, we observe the dynamics for a mild grooming behaviour where αL = 0.4,αN = 0.6,αA = 0.5 and κ = 0.1 × 10−5. The equilibrium value for all stages are higher when there is no grooming behaviour. In particular, the value of the adult egg laying females at the equilibrium is 693 for a mild resistance behaviour and 1.9×103, when there is no resistance (Figure 3 and the left side of Figure 6). We also see that by decreasing the resistance solutions 78 M. ALAVINEJAD, J. SADIKU, AND J. WU (a) (b) (c) (d) Figure 2. Case 1, Rv0 < 1 where βL = 0.6×10−4, βN = 1.8×10−4 and p = 200 yields Rv0 = 0.89. with non-oscillatory behaviour show damped oscillation. In a maximum intensified grooming behaviour the tick attachment rates to hosts with resistance are reduced to 0, therefore high resistance of hosts affects the tick equilibrium values significantly. For instance, in Figure 6 the equilibrium value for Aelf reduces from 1.9 × 103, when there is no resistance, to 78 when the resistance is very high. Comparing the right side of Figure 5 with 7, demonstrates the effect of resistance factors on the dynamical behaviour of the solutions. Reducing the resistance from high to a mild resistance results in an increase in the value of the equilibrium of Aelf from 78 to 1600. However, in the absence of host resistance, the tick population at different stages oscillate about a positive equilibrium (Aelf ≈ 2.7×103). In other words, by decreasing the grooming behaviour (increasing the value of αL, αN and αA from 0 to 1), there is more available resources for ticks to feed on. Therefore, the dynamical behaviour of tick population at different stages changes from solutions converging to the positive equilibrium to oscillatory solutions. The dynamics of the feeding ticks are similar to those of questing ticks and therefore we exclude the pictures on this paper. When we ignore the resistance behaviour in case 2 and 3, the host population with resistance Hr+ is equal to 0 and it reaches a positive equilibrium point when αL = αN = αA = 0. 4.2. LHS and PRCC. We perform Latin Hypercube Sampling to further analyze the effects of each parameter on the dynamics of each life stage of the ticks [1, 8] Before we proceed to performing PRCC a verification of monotonicity is necessary to ensure the correct range of the parameters for PRCC analysis. Next, we calculate the PRCC, which determines the contribution of each parameter to the output variable such the population of larvae questing. A PRCC value significantly greater than 0 indicates a positive correlation and for PRCC significantly less than 0, a negative correlation between the parameter and the output [14]. In Figure 8, the PRCC for the larvae questing population demonstrates MODELING THE IMPACT OF HOST RESISTANCE ON STRUCTURED TICK POPULATION DYNAMICS 79 (a) (b) (c) (d) Figure 3. In Case 2 the values of p and κ have changed to p = 1500 and κ = 0.1×10−5 and the reproduction number increased to Rv0 = 6.71. The simulations run for a time span of 10000 days. The equilibrium points for each stage of questing, feeding and adult egg laying female tick are as follows: Lq = 6.5 × 107,Nq = 1.6 × 106,Aq = 1.6 × 105 Lf = 2.9×106,Nf = 2.9×105,Af = 1.4×105,Aelf = 693. In addition, the equilibrium point of the host with resistance is 13. Figure 4. The solutions oscillate about the equilibrium point as we change the value of αA in the interval [0, 1] for parameter values in case 3. The values for αL and αN are 0.6 and 0.8. 80 M. ALAVINEJAD, J. SADIKU, AND J. WU (a) (b) (c) (d) Figure 5. In Case 3 the values of βL and βN have changed to βL = 1.2×10−4,βN = 3 × 10−4 producing a higher reproduction number, Rv0 = 16.9. The simulation are again running for a time span 10000 days. The equilibrium points for each stage of questing, feeding and adult egg laying female tick are as follows: Lq = 5.7 × 107,Nq = 2.3 × 106,Aq = 3.8 × 105,Lf = 4.8 × 106,Nf = 6.9 × 105,Af = 3.2 × 105,Aelf = 1600. In addition, the equilibrium point of the host with resistance is 14. the negative correlation with the death rates dAelf , dNm, dLm, dNq, dAq, dLq, dE and dLq having the highest effect on this stage. The detachment rate δ does not have a an impact, however the parameters related ticks’ biological characteristics, p, q, �, have a significant effect. We also observe that the host finding rates βA, βN ,βL, have positive correlation with larvae questing dynamics. For the values of most parameters that are taken from the literature, we would expect to see a reasonable correlation between the parameter and the output (in a range where the output is monotonically increasing or decreasing with parameter). For instance the output value of Lq (and therefore Lf ) at the equilibrium is supposed to decrease with an increase of the larvae questing death rate (negative correlation). 5. Conclusion In this paper we formulated a delay differential model for black leg ticks, stratified based on stage and activity, with a particular focus on the host grooming behaviour. The basic reproduction number was calculated and the condition for local stability of tick-free equilibrium, for which the tick population go extinct, and also for existence and uniqueness of a positive equilibrium was given. Model param- eterization and numerical simulations were carried out to demonstrate the dynamics of tick and host population with and without the grooming behaviour and the effect of the resistance factor on the value of equilibrium points are studied. Parameters related to the grooming and resistance factors, αL, αN , αA, and κ have no effect on the initial growth rate of ticks since these parameters do not change the MODELING THE IMPACT OF HOST RESISTANCE ON STRUCTURED TICK POPULATION DYNAMICS 81 (a) (b) (c) (d) Figure 6. The parameter values are the same as in Case 2 except the αL = αN = αA = 1 (on the left). The equilibrium points are as follows: Lq = 5.0 × 107, Nq = 2.3×106, Aq = 2.4×105, Aelf = 1.9×103. There is no resistance and hence Hr+ = 0. In case of αL = αN = αA = 0 (on the right) the equilibrium points are Lq = 1.3×107, Nq = 3.2 × 105, Aq = 2.4 × 104, Aelf = 78. Since now we introduce resistance, Hr+ = 10. value of Rv0 . However, with an increase of the intensity of the grooming behaviour from no resistance to a high level of resistance, where either the hosts show intensified grooming behaviour or ticks are withdrawn from feeding or dead, the values of equilibrium points of all tick stages decrease. From the numerical simulations we observed structural changes of the dynamical behaviour of the tick popula- tion by changing the parameter values reflecting the effect of the host resistance. Also, the intensified resistance results in higher equilibrium values for Hr+. A sensitivity analysis of the positive equilibrium value to the parameters was carried out by perform- ing LHS and PRCC. From PRCC we observed high positive correlation between the maximum number of eggs per female adult tick (p) and larvae questing; as more eggs are produced the higher the number of larvae questing. The female proportion parameter (�) is also positively correlated to larvae questing. As the female rate proportion increases the higher number of egg production and therefore increasing the value of larvae questing. In contrast, the value of the strength of density dependence (q) and death rate of larvae questing (dLq) are negatively correlated with the population of larvae questing. As the death rate increases there will be a lower population size of larvae questing. Lastly, as the number of larvae questing increases there will be harder to find resources to survive, hence as q increases the number the Lq decreases. This study has some limitations. The death rates are assumed to be constants for each stage of the tick and we have ignored the possibility of death during the feeding process resulting from serous exudes 82 M. ALAVINEJAD, J. SADIKU, AND J. WU (a) (b) (c) (d) Figure 7. The parameter values are the same as in Case 3 except αL = αN = αA = 1 (the left). The equilibrium points are as follows: Lq ≈ 2.8 × 107, Nq ≈ 2.1 × 106, Aq ≈ 3.5×105, Aelf ≈ 2.7×103. Since resistance factor is not introduced the Hr+ = 0. On the right side the αL = αN = αA = 0 and the equilibrium points are as follows: Lq = 1.4 × 107, Nq = 4.0 × 105, Aq = 3.8 × 104, Aelf = 80. The resistance factor increase the population size from zero to Hr+ = 11. Figure 8. PRCC for most of the parameters used in the model at the equilibrium point of Lq. The value of each parameter is taken from 1and Case 2 for a range of (+/−)20% MODELING THE IMPACT OF HOST RESISTANCE ON STRUCTURED TICK POPULATION DYNAMICS 83 which could engulf the tick. Also, interpreting the host resistance as a kind of immunity to ticks we can consider the situation where the host resistance decreases in time the hosts lose immunity to ticks. The molting process is demonstrated by constant delay functions. Future work could incorporate the temperature and humidity on molting process and explore further the effects on tick dynamics. References [1] S. M. Blower and H. Dowlatabadi, Sensitivity and uncertainty analysis of complex models of disease transmission: an HIV model, as an example, International Statistical Review/Revue Internationale de Statistique 92(1994), 229-243. [2] J. Bowessidjaou, Effects and duration of resistance acquired by rabbits on feeding and egg laying in Ixodes Ricinius, L. Experientia 33(1977), 528-530. [3] S. J. Brown, Highlights of contemporary research on host immune response to ticks, Elsevier Science Publishers 28(1988), 321-334. [4] S. J. Brown and P. W. Askenase, Rejection of ticks from guinea pigs by anti-hapten-antibody-mediated degranulation of basophils at cutaneous basophil hypersensitivity sites: role of mediators other than histamine, The Journal of Immunology 134(1985), 1160-1165. [5] [S. J Dumler, Molecular diagnosis of Lyme disease: review and meta-analysis, Molecular Diagnosis 6(2001), 1-11 [6] G. Fan, H. R Thieme, and H. Zhu, Delay differential systems for tick population dynamics, Journal of Mathematical Biology, 71(2015), 1017-1048. [7] H. D. Gaff and L. J. Gross, Modeling tick-borne disease: A metapopulation model, Bulletin of Mathematical Biology 69 (2007), 265-288. [8] B. Gomero, Latin Hypercube Sampling and Partial Rank Correlation Coefficient Analysis Applied to an Optimal Control Problem, MSc Thesis, University of Tennessee at Knoxville, 2012. [9] R. Jennings, Y. Kuang, H. R. Thieme, J. Wu, and X. Wu, How ticks keep ticking in the adversity of host immune reactions, Journal of Mathematical Biology 78(2019), 1331-1364. [10] R. C. Johnson, G. P. Schmid, F. W. Hyde, A. G. Steigerwalt and D. J. Brenner, Borrelia burgdorferi sp. nov.: Etiologic agent of Lyme disease, International Journal of Systematic and Evolutionary Microbiology 34(1984), 496- 497. [11] Y. Lou and J. Wu, Tick seeking assumptions and their implications for lyme disease predictions, Ecological Com- plexity 17(2014), 99-106. [12] W. Lumsden, Advances in parasitology, Vol. 18, Academic Press, 1980. [13] N. K. Madhav, J. S. Brownstein, J. I. Tsao and D. Fish, A dispersal model for the range expansion of blacklegged tick (acari: Ixodidae), Journal of Medical Entomology 41(2004):842-852. [14] S. Marino, I. B. Hogue, C. J. Ray and D. E. Kirschner, A methodology for performing global uncertainty and sensitivity analysis in systems biology, Journal of Theoretical Biology 254(2008), 178-196. [15] N. H. Ogden, M. Bigras-Poulin, C. J. O’callaghan, I. K. Barker, L. R. Lindsay, A. Maarouf, K. E. Smoyer-Tomic, D. Waltner-Toews and D. Charron, A dynamic population model to investigate effects of climate on geographic range and seasonality of the tick Ixodes Scapularis, International Journal for Parasitology 35(2005), 375-389. [16] J. Piesman, T. N. Mather, R. Sinsky and A. Spielman, Duration of tick attachment and borrelia burgdorferi trans- mission, Journal of Clinical Microbiology 25(1987), 557-558. [17] Public Health Agency of Canada: Surveillance of Lyme Disease, https://www.canada.ca/en/public- health/services/diseases/lyme-disease/surveillance-lyme-disease.html [18] S. Randolph, Tick ecology: processes and patterns behind the epidemiological risk posed by Ixodid ticks as vectors, Parasitology 129(S1) (2004), S37-S65. [19] R. Rosa and A. Pugliese, Effects of tick population dynamics and host densities on the persistence of tick-borne infections, Mathematical Biosciences 208(2007), 216-240. [20] R. Rosa, A. Pugliese, R. Norman and P. J. Hudson, Thresholds for disease persistence in models for tick-borne infections including non-viraemic transmission, extended feeding and tick aggregation, Journal of Theoretical Biology 224(2003), 359-376. [21] H. L Smith, Monotone semifows generated by functional differential equations, Journal of Differential Equations 66(1987), 420-442. [22] A. C. Steere, J. Coburn and L. Glickstein, The emergence of lyme disease, The Journal of Clinical Investigation 113(2004), 1093-1101. [23] A. C Steere, R. L. Grodzicki, A. L. Kornblatt, J. E. Craft, A. G. Barbour, W. Burgdorfer, G. P. Schmid, E. Johnson and S. E. Malawista, The spirochetal etiology of Lyme disease, New England Journal of Medicine 308(1983), 733-740. 84 M. ALAVINEJAD, J. SADIKU, AND J. WU [24] B. M. Wagland, Host resistance to cattle tick (boophilus microplus) in brahman (bos indicus) cattle II. The dynamics of resistance in previously unexposed and exposed cattle, Australian Journal of Agricultural Research 29 (1978), 395- 400. [25] X. Wang and X. Zhao, Dynamics of a time-delayed lyme disease model with seasonality, SIAM Journal on Applied Dynamical Systems 16(2017), 853-881 [26] X. Wu, Modeling the impact of climate change on tick population dynamics and tick-borne disease spread, PhD Thesis, York University, 2013. [27] X. Wu, F. M. Magpantay, J. Wu and X. Zou, Stage-structured population systems with temporally periodic delay, Mathematical Methods in the Applied Sciences 38(2015), 3464-348. [28] X. Wu, V. R. S. K. Duvvuri and J. Wu, Modelling dynamical temperature influence on tick Ixodes Scapularis population, International Congress on Environmental Modelling and Software 529, 2010. Corresponding author. Department of Mathematics and Statistics, York University E-mail address: mahnazal@yorku.ca Department of Mathematics and Statistics, York University E-mail address: jemisa18@yorku.ca Department of Mathematics and Statistics, York University E-mail address: wujh@yorku.ca