MODELING THE IMPACT OF HOST RESISTANCE ON STRUCTURED TICK POPULATION DYNAMICS

. 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 reﬂecting the full biological stages of ticks. Here we develop and simulate a delay diﬀerential equation model, with a particular focus on resistance resulting in grooming behaviour. We calculate the basic reproduction number using the spectral analysis of delay diﬀerential 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 diﬀerent stages.


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 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 ticksensitized 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 transmission 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 persistence 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.

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 A elf , adult questing (A q ) and feeding (A f ). 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: H r+ represents the bitten host compartment that have developed with immunity; H r− 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 densitydependent 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. Per capita mortality rate of L q 0.6 × 10 −2 per day [15] d Lm Per capita mortality rate of L m 0.3 × 10 −2 per day [15] d Nq Per capita mortality rate of N q 0.6 × 10 −2 per day [15] d Nm Per capita mortality rate of N m 0.2 × 10 −2 per day [15] d Aq Per capita mortality rate of A q 0.6 × 10 −2 per day [15] d A elf Per capita mortality rate of A elf 1 per day [28] d E 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  [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   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: where γ(A) = pAe −qA is the birth function. Here, we use the following equation for the host population dynamics where H(t) = H r− (t) + H r+ (t). Note that the positive equilibrium of this equation is given byH = Interpreting K as an environmental constraint, and in order to haveH ≤ K we assume c ≥ (b − µ), withH = K when the equality holds.
Let C τ be the product of Banach spaces .
Then the right hand side of the Equation (2.1) is given by (2.5) We assume the initial data is non-negative. So we will assume the initial data X 0 is in the Banach Also, for the initial data to be continuous and positive we assume: (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.
Since initial data for A elf on [−τ, 0] is positive, the derivative of L q at t 1 is positive and therefore L q (t) is increasing, so it can not be negative. The same argument can be applied for [τ, 2τ ]. This proves that L q (t) ≥ 0 for all t ≥ 0. If there exists t 2 such that L f (t 2 ) = 0, then the derivative of L f at t 2 is given by which is positive since L q (t 2 ), H r+ (t 2 ) and H r− (t) are positive. Thus L f is increasing at t 2 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 H r− and H r+ 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) Using the fact that γ(A elf (t)) ≤ p/qe for all t ≥ 0, we see that Integrating the next equations and taking the supermom we have: 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 atH = (b − µ)K/c, the limiting system is as follows From now on we refer to this system as the main system of our model unless otherwise stated.

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.
3.1. Equilibria. Let X * denote the vector (L q , L f , N q , N f , A q , A f , A elf , H r+ ) in R 8 , 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: At the tick-free equilibrium, where all tick stages are equal to zero, we have We want to derive conditions for existence and uniqueness of a (strongly) positive equilibrium point (X i > 0 for all i = 1, · · · , n). From the equations in (3.1) we get the following N ) and s 3 = ψe −d Nm τ (N,A) . From the first equation in the system (3.1) we get and therefore Since γ(A elf ) = pA elf e −qA elf we have the following cases: A elf = 0 or (3.5) Finally, we reduce the system (3.1) to the following system given that H r+ =H and Ω(H r+ ) = 0 (it can be proved that this holds). Substituting this in the equation  (i) Γ is a rational function and is strictly increasing for 0 < H r+ <H; (ii) Γ(0) > 0 is given by (iii) G is a negative exponential function and it approaches zero (exponentially) as H r+ approaches H; (iv) G(0) = p.
From these properties we can see that the equation (3.7) has at least one solution 0 < H r+ <H, if and only if G(0) > Γ(0), i.e., If R v 0 > 1, then system (2.9) has a positive equilibrium point. If additionally d dH r+ holds, then the positive equilibrium is unique.

3.2.
Stability of the tick-free Equilibrium. First we linearize System (2.9) about a given equilibrium point using the Fréchet derivative of the function F (X), given by the right hand side of the System (2.9): The linearized system is given by The linearized system about the tick-free equilibrium point is as follows: 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. Theorem 3.2. If R v 0 < 1, then X = 0 is the only equilibrium point of the system (2.9) and is locally asymptotically stable. When R v 0 > 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): The matrix J can be written as J = F − V . The zero equilibrium is locally asymptotically stable if ρ(F V −1 ) < 1 (ρ is the spectral radius of F V −1 ) and it is unstable if ρ(F V −1 ) > 1. We can see that

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., R v 0 < 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 × 10 3 , when there is no resistance (Figure 3 and the left side of Figure 6). We also see that by decreasing the resistance solutions 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 A elf reduces from 1.9 × 10 3 , 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 A elf from 78 to 1600. However, in the absence of host resistance, the tick population at different stages oscillate about a positive equilibrium (A elf ≈ 2.7 × 10 3 ). 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 H r+ is equal to 0 and it reaches a positive equilibrium point when α L = α N = α A = 0.

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 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 R v 0 = 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: L q = 6.5 × 10 7 , N q = 1.6 × 10 6 , A q = 1.6 × 10 5 L f = 2.9×10 6 , N f = 2.9×10 5 , A f = 1.4×10 5 , A elf = 693. In addition, the equilibrium point of the host with resistance is 13.  the negative correlation with the death rates d Aelf , d N m , d Lm , d N q , d Aq , d Lq , d E and d Lq 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 L q (and therefore L f ) at the equilibrium is supposed to decrease with an increase of the larvae questing death rate (negative correlation).

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 parameterization 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 (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: L q = 5.0 × 10 7 , N q = 2.3 × 10 6 , A q = 2.4 × 10 5 , A elf = 1.9 × 10 3 . There is no resistance and hence H r + = 0. In case of α L = α N = α A = 0 (on the right) the equilibrium points are L q = 1.3 × 10 7 , N q = 3.2 × 10 5 , A q = 2.4 × 10 4 , A elf = 78. Since now we introduce resistance, H r+ = 10.
value of R v 0 . 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 population by changing the parameter values reflecting the effect of the host resistance. Also, the intensified resistance results in higher equilibrium values for H r+ .
A sensitivity analysis of the positive equilibrium value to the parameters was carried out by performing 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 (d Lq ) 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 L q 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  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.