A coupled algebraic-delay differential system modeling tick-host interactive behavioural dynamics and multi-stability

We propose a coupled system of delay-algebraic equations to describe tick attaching and host grooming behaviors in the tick-host interface, and use the model to understand how this tick-host interaction impacts the tick population dynamics. We consider two critical state variables, the loads of feeding ticks on host and the engorged ticks on the ground for ticks in a particular development stage (nymphal stage) and show that the model as a coupled system of delay differential equation and an algebraic (integral) equation may have rich structures of equilibrium states, leading to multi-stability. We perform asymptotic analyses and use the implicit function theorem to characterize the stability of these equilibrium states, and show that bi-stability and quadri-stability occur naturally in several combinations of tick attaching and host grooming behaviours. In particular, we show that in the case when host grooming is triggered by the tick biting, the system will have three stable equilibrium states including the extinction state, and two unstable equilibrium states. In addition, the two nontrivial stable equilibrium states correspond to a low attachment rate and a large number of feeding ticks, and a high attachment rate and a small number of feeding ticks, respectively.


Introduction
Tick-borne diseases (TBDs) have been imposing significant public health challenges globally. The TBD transmission relies on the tick-host interaction where ticks may acquire and/or transmit the pathogen from the host by taking a blood meal from the host. There are systemic transmission and co-feeding routes, both depending on the cooccurrence of ticks (specially ticks at two different stages for co-feeding transmission to take place) on the same host, so understanding the tick load distribution dynamics over the host and its implication for the tick population dynamics, the main focus of this study, is important. The tick load distribution process over the host is a dynamical process, governed by the tick attachment/fixation and host grooming behaviours. Our study shows that these tick and host individual behaviors, the host response to tick attaching and/or to tick feeding, and their combinations can yield a complicated tick-host interaction leading to multi-stability where tick densities can converge to a tick extinction, or a lower level tick persistence, or a higher level tick persistence equilibrium state depending on the initial conditions.
Tick life cycle includes four stages: egg, larvae, nymph and adult. Larval and nymphal ticks seek blood meals from small rodents, like mice and bird, to molt into nymphal and adult stage, respectively. While it is hard to find larval and nymphal ticks on the hosts since they are small [only around 1-2 mm in size (Lindquist and Vapalahti 2008)], ticks at these stages are very important for the tick-borne pathogen transmission as vertical transmission (from egg-laying infected ticks to eggs) are limited. Adult ticks prefer questing large mammals, such as deer and domestic livestock, both infected nymphal and adult ticks can bite on human, passing the pathogen to the human. After the final blood meal, adult tick will lay egg to complete the life cycle from eggs to egg-laying ticks.
Important to the survival of an egg through the life cycle is the attachment/fixation success of the questing tick. Existing models have assumed a constant attachment success rate, so the fixation rate of an engorged tick to a host for feeding to develop into the next stage is the questing rate times the attachment/fixation success rate. However, it was reported (Voordouw 2015;Wang et al. 1998Wang et al. , 2001Rechav and Nuttall 2000) that ticks may pool their saliva to enhance the immunomodulatory manipulation of the host organism. The resulted cooperative feeding could increase the cost-benefit ratio of resource extraction from the host relative to per capita investment in tick saliva production. Reflecting this cooperative feeding requires a (feeding tick) densitydependent attachment/fixation rate, as we will consider in our model formulation.
Equally important to the survival of an egg through the life cycle is the grooming behavior of the host, leading to (feeding tick) density-dependent grooming rate (or feeding tick survival rate) in some of the proposed mathematical models. Pioneering observation of the phenomenon called acquired tick resistance (ATR), by Trager (1939), showed that upon repeated tick infestations, hosts develop an immune response to derail subsequent tick challenges, and tick-immune hosts rapidly reject ticks within the first 24 h of tick attachment. See recent reviews (Narasimhan et al. 2021;Yoshikawa et al. 2020) for the current knowledge of ATR and key events in the tick-host interaction to enable or disable tick feeding. Evidence provided by Hart and his colleagues (Hart 2000) provided support to "the concept that the delivery of bouts of grooming reflects programmed grooming", namely, "grooming occurs in response to an endogenous generator that produces grooming bouts at periodic intervals, resulting in removing of ticks before they attach and begin to feed". Correspondingly, in our model formulation, we will consider (feeding tick) density-dependent grooming and/or attachment/fixation-dependent grooming rate.
In summary, a mathematical description of the tick population dynamics must take consideration of the tick attachment and fixation dynamics, and host grooming dynamics. Here, we develop a coupled system involving two state variables, the number of engorged nymphal ticks and the number of feeding nymphal ticks. We start with an evolution equation for the feeding tick density (with respect to the time since feeding nymphal ticks attach to the host) governed by the density-dependent grooming rate, and subject to an initial condition involving density-dependent attachment and fixation rate. Integration along a characteristic equation leads to an algebraic equation for the total feeding (nymphal) ticks. We then couple this algebraic equation with a delay differential equation for the engorged nymphal tick dynamics. Since questing nymphal ticks come from the engorged nymphal ticks with a delay after further development and production, the coupled system becomes a closed feedback system with delay.
The model will be formulated in Sect. 2; the model's equilibrium structure is described in Sect. 3; and the stability of equilibrium states is discussed in Sect. 4 using a perturbation argument since the feeding duration is relatively short in comparison with the life cycle. Additional discussions on how this tick population dynamics can be extended to model tick-borne disease transmission dynamics involving the co-feeding transmission route are provided in the final section.

The algebraic-delay coupled system
We focus on the developmental stage, nymphal stage, where questing and attaching rate of the ticks and the grooming rate of the hosts may depend on the feeding nymphal tick loads on the host.

The model derivation and simplification
Let Q(t), F(t) and E(t) be the numbers of questing nymphal ticks, feeding nymphal ticks and engorged nymphal ticks, respectively, at time t. We make it explicit our standing assumptions: (i) The attachment rate of questing nymphal ticks is a function ρ(F(t)) of the number of feeding nymphal ticks F(t). This is to reflect the fact that the attachment success of questing ticks depends on the total amount of feeding nymphal ticks on the hosts. As the total number of hosts for nymphal ticks is relatively static, the average feeding ticks on the hosts is proportional to F(t). In what follows, we assume ρ : [0, ∞) → [0, ∞) is a C 1 -smooth function. Moreover, ρ(F) > 0 if F > 0. (ii) The drop off rate of feeding nymphal ticks is a function ν(F(t)) of the number of feeding nymphal ticks F(t). This is to describe the host grooming behaviors, and can also be used to describe the cooperative co-feeding behaviours of the nymphal ticks (that enhances the immune-induced feeding). We also assume that Let n(t, a) be the density of feeding nymphal ticks at time t with respect to feeding duration a since they attach to the host. Then we have the structured feeding tick population dynamics model where T is the average feeding duration of nymphal ticks. In this formulation, the evolution describes the grooming dynamics while the boundary condition describes the attaching behaviours. The structured population model can be easily solved using the method of integration along characteristics (Evans 1998). Namely, let then we can rewrite (1) into the following form Integration of (3) yields which, from the definition of (2), is equivalent to Setting t − s = a, we obtain Clearly, the total amount of feeding nymphal ticks can be expressed as the integral of n(t, a) over the feeding interval, i.e., which is an implicit equation for F(t). Then the dynamics for engorged nymphal ticks follows where δ is the exit rate of engorged nymphal ticks. We now link the questing nymphal ticks at the current time to engorged nymphal ticks in the past through the life cycle. Let η 1 be the survival probability from engorged nymphal ticks to adult egg-production ticks, σ be the egg production rate, and η 2 be the survival probability from eggs to questing nymphal ticks. Assume that τ 2 represents the delay from eggs to questing nymphal ticks through the necessary developments and τ 1 is the delay from engorged nymphal ticks to adult egg-laying ticks. Hence, questing nymphal ticks at time t take the following form For simplicity, we denote τ = τ 1 +τ 2 and η = η 2 σ η 1 . From the definition of feeding duration T of nymphal ticks, it is easy to see that the life cycle of tick population is T + τ . We obtain the following coupled system to describe the dynamics of feeding and engorged nymphal ticks (4)

Fundamental theory
We now show that the above coupled system of differential-algebraic system is equivalent to a coupled system of delay differential equations subject to a matching condition of the initial data. Differentiating the right-hand side of the algebraic equation for F(t), we get This is equivalent to the algebraic equation if the following matching condition is met: Now, we can use the fundamental theory for functional differential equations (Hale 1977) to conclude that for any , there is one and only one solution of the following coupled system of delay differential equations The solution (E(t), F(t)) ∈ R 2 is defined for t ≥ 0. With the matching condition (6), we conclude from (5) that and hence F satisfies the algebraic equation in system (4). We now show that if φ(θ), ψ(θ) > 0 for θ ∈ [−τ − T , 0], and if the matching condition (6) hold, then E(t), F(t) ≥ 0 for all t ≥ 0. To prove this, we use the continuous dependence of solutions on parameter > 0 for the following system: Denote the solution of (7) by ( not true for all t ≥ 0, then there must be the first t * ≥ 0 such that E (t * ) = 0 and E (t) > 0, F (t) > 0 for all t ∈ [0, t * ). Therefore, we have d dt E (t)| t=t * ≤ 0. But using the first equation of (7), we yield That is a contradiction. Thus, we have the following existence-uniqueness, and positiveness result.
, and if the matching condition (6) is satisfied, then system (4) has one and only one solution defined for all t ≥ 0. This solution is non-negative, namely, Remark 1 There are two other approaches to establish the fundamental theory for the well-posedness of the coupled system we formulated. First of all, we can solve the algebraic equation by using the implicit function theory to obtain t] ) and then substitute this to the first equation to obtain a single functional differential equation for E(t) although the right hand side is given implicitly. Alternatively, we can rewrite the algebraic equation as Therefore, the coupled system can be regarded as a special case of the neutral functional differential equation , the neutral operator D : X → R 2 and the functional f : X → R 2 are given by The fundamental theory of neutral functional differential equations including the principle of linearization can be found in Hale (1977), Hale and Lunel (1993). See also Barbarossa et al. (2014) and Gourley and Kuang (2009) for neutral equations arising from structured population dynamics in other settings.

Equilibria and stability of trivial state
The equilibrium of model (4) satisfies the following nonlinear equations

Trivial equilibrium and its stability
Clearly, model (4) always has a trivial equilibrium P 0 (0, 0). We linearize model (4) at the zero equilibrium P 0 to obtain The stability of the trivial equilibrium P 0 is determined by the following characteristic equation Therefore, the trivial equilibrium P 0 is locally asymptotically stable if δ > ρ(0)ηe −ν(0)T , and unstable if δ < ρ(0)ηe −ν(0)T . This is because of the positive feedback in model (9), and the semi-group generated by this linear delay differential equation is order-preserving, and the stability of the zero solution is the same as that of the following ordinary differential equation by using the monotone dynamical system theory (Smith 1987(Smith , 1995. The linearization of the coupled system for the algebraic equation is Therefore, F(t) → 0 as t → ∞.

Non-trivial equilibria
The nontrivial equilibrium where δ = ηρ(F * )e −ν(F * )T < ηρ(F * ). Therefore, the existence and multiplicity of nontrivial equilibrium depends on behaviors of tick attachment rate ρ(F) and the host grooming rate ν(F). We consider several scenarios of the tick-host interface.

Constant attaching and grooming
We first consider the simplest case of constant attachment rate and grooming rate, with ρ(F) = p and ν(F) = μ 0 , where p, μ 0 > 0 are positive constants. Model (4) becomes a linear system as no nonlinearity involves: The first equation is a scalar delayed differential equation with a positive delayed feedback. Clearly, the basic reproduction number is R 0 = pηe −μ 0 T δ −1 , obtained from the multiplication of reproduction and survival probability during the life cycle except the nymphal tick engorgement with sojourn time δ −1 . An application of the Krein-Rutman theorem (see Smith 1987Smith , 1995 shows that solution of E(t) with a non-trivial non-negative initial value on [−τ −T , 0] is convergent to 0 or ∞ as t → ∞ when R 0 < 1 and R 0 > 1 respectively. Correspondingly, using the second (integral) equation, we obtain that Remark 2 Note also that (10) has infinitely many positive equilibria in the critical case when R 0 = 1. Namely, when there are infinitely many positive equilibria (11). Then the first differential equation of (11) is reduced into This type of delay differential equation was studied previously in Haddock-Terjeki (1983) and it has the so-called asymptotic constancy property. That is, lim t→∞ E(t) = E c exists and is a constant for each given solution. To determine the value E c for each given solution, we first observe that from which it follows that Therefore

Density-dependent monotone attaching and grooming rates
A more biologically realistic situation is when the attachment rate decreases and grooming rate increases with tick loads on the host. We consider the prototypical case where c > 0, μ > 0 and μ 0 > 0 are constants. In this case, we can reduce model (4) into Model (13) has a unique positive equilibrium if R 0 > 1, and no positive equilibrium to conclude that E(t) → 0 as t → ∞, and then use the integral inequality

Cooperative feeding and density-dependent grooming
Recall that ticks may pool their saliva to enhance the immunomodulatory manipulation of the host organism. The resulted cooperative feeding could increase the cost-benefit ratio of resource extraction from the host relative to per capita investment in tick saliva production. Therefore, in cooperative feeding, the attachment/fixation rate is an increasing function of the feeding tick density on the host. However, since ticks prefer seeking for soft and thin areas of host skin that are well-supplied with blood, there is a maximum capacity to accommodate tick attachments. To describe this cooperative and self-limiting feeding attaching process, we consider the case where the attachment/fixation rate ρ(x) is an initially increasing function that becomes decreasing after the capacity (c) is reached: there exist two constants p > 0 and c > 0 satisfying ρ (0) = p > 0, ρ (F) > 0 for F ∈ (0, c) and ρ (F) < 0 for F ∈ (c, +∞). For the sake of simplicity, we use the Ricker function (see Ricker 1975) as a prototypical attachment rate function, i.e., We will couple this cooperative and self-limiting attachment with the densitydependent grooming rate: responding to increasing feeding ticks, the host groom more frequently, resulting in dropping off rate function ν(F) being an increasing function, i.e., Then model (4) can be rewritten as (15) To look at a positive equilibrium F of model (15), we consider positive solutions of the first equation of (10), namely, The function g 1 (x) := xe −( 1 c +μT )x changes its monotonicity (from increasing to decreasing once), and researches its maximum where F * ± are two positive solutions of (16) and the corresponding E * ± are given by the second expression of (10) with F being replaced by F * ± respectively.

Grooming reactive to tick biting
When hosts groom in response to tick biting, we have with f (ρ) being an increasing function of ρ ≥ 0. Here, we consider a linear function f (ρ) = kρ. Then model (4) becomes The first equation of (10) for the equilibrium state becomes an equation for ρ: Let g 2 (ρ) := ρe −kT ρ . Then g 2 (ρ) = e −kT ρ (1 − kT ρ), and we obtain the maximum value (ekT ) −1 of the function g 2 (ρ) at ρ = (kT ) −1 . For a corresponding attachment rate function, we use the function where r > 0 and j > 0. The choice of this function is motivated by the Holling type III functional response to describe the improvement of the tick's searching for an appropriate soft and thin area of the host skin. Note that we have separated the ticks on the hosts into two classes: those questing ticks who have successfully attached to the host but not yet feeding on the host, and those who are feeding, therefore the questing functional response function is ρ(F)Q, which would be the saturation function r F 2 /(1 + j −2 F 2 ) should Q = F. We use r for the "attack" rate of the type III functional response (normally denoted by a, but this has been reserved for the age-variable in our current study). The parameter j is related to the attach rate (r ) and the handling time (normally denoted by h), that is, j −1 = rh. The function ρ(F) reaches its maximum value r j/2 at F = j. A simple geometric argument yields Theorem 3 For the combination of the attachment rate (18) and the linear grooming behaviour f (ρ) = kρ, we have the following results about the equilibrium structure: (i) When δη −1 > (ekT ) −1 , there is no positive equilibrium in model (17); (ii) When δη −1 < (ekT ) −1 , the equation g 2 (ρ) = δη −1 has two positive solutions ρ * − and ρ * + , which satisfy that ρ * − < (kT ) −1 < ρ * + . We have the following subcases: are two positive solutions of the equation ρ(x) = ρ * − ; and the corresponding E * −,± are given by the second expression of (10) with F being replaced by F * −± respectively; Fig. 1) are positive solutions of the equations ρ(x) = ρ * − and ρ(x) = ρ * + , respectively, and the corresponding E * ±,± are given by the second expression of (10) with F being replaced by F * ±± respectively.
We illustrate the case (ii3) in Fig. 1, where we see the model has four positive equilibria.

Stability and characteristic equation
We translate the positive equilibrium P + into the origin by the translation Model (19) is a delayed differential-algebraic system, also called a degenerate differential system with delay. Linearization of system (19) takes the following form Looking for a non-trivial exponential vector function as a solution x(t) = e λt x 0 , where (x 0 , y 0 ) T = 0, we obtain where In order for (x 0 , y 0 ) T to be a nonzero vector, the determinant of the coefficient matrix of the system (22) should be zero, i.e., the characteristic equation of (20) takes the following form It is highly nontrivial to describe the distribution of zeros of this characteristic equation in the general case. In what follows, we take advantage of the fact that the feeding duration T is relatively small (a few days) comparing with the tick life cycle in the natural world that can be as long as 3 years.
We start with the simplest case where = 0. If c/e > δ/(ηp), this special case reduces the positive equilibria of model (15) are the solutions of the following equation from which we obtain the corresponding E-coordinate: Using the implicit function theorem, when T = > 0, it can be seen that model (15) has two positive equilibria (E * − ( ), F * − ( )) and (E * + ( ), F * + ( )), where F * − ( ), F * + ( ) satisfying F * − ( ) < c < F * + ( ) are the two positive solutions of the following equation Now, considering T = as a variable parameter, we focus on searching for the characteristic eigenvalue of linearized system (20) with respect to . As the equilibrium E( ) has a singularity as → 0, the characteristic equation has a singularity at = 0 so we need to single out this singularity. Define Therefore, we consider G(λ, ) = 0 as the characteristic equation of system (19). In what follows, we will calculate the series expansion of G(λ, ) = 0 with respect to . It is clear to see that Thus, a 12 (λ, ) → ηδ −1 ρ (F 0 ± )F 0 ± and a 22 (λ, ) → 0, as → 0.

Suppose that the series expansions of F( ) and E( ) with the following forms
where o( ) represents higher order terms. From the expressions of ρ(F) and ν(F), we obtain the following series expansions Combining the above series expansions, the fist equation of (10) can be rewritten as which leads to δ = ηρ 0 and ρ 1 = ν 0 ρ 0 .
Based on the relationship in (28), we have Using the second equation of (8), we have Substituting the series expansions of F( ), E( ), ρ(F), ν(F) into (29), we yield (27) and the relationship in (28), we have

From all coefficients in
Thus, the series expansions in (26) can be expressed into Assume that the series expansion of eigenvalue satisfying G(λ, ) = 0 is given by Based on the expressions of a i j in (22), we can calculate their series expansions as follows The characteristic polynomial G(λ, ) becomes From the characteristic equation G(λ, ) = 0, it follows that Using the expressions in (27) and (28), we can simplify (30) into the following form The stability of the two positive equilibria (F * − ( ), E * − ( )) and (F * + ( ), E * + ( )) depends on the sign of the first non-zero term λ 0 . In the series expansion of F( ) in (26), we know that the first term F 0 is equal to the equilibrium F * − ( ) or F * + ( ). Since F * − ( ) < c < F * + ( ), it follows that Sign(λ 0 ) > 0, at the positive equilibrium (F * − ( ), E * − ( )), Sign(λ 0 ) < 0, at the positive equilibrium (F * + ( ), E * + ( )).
Then we have the following result: Theorem 4 Assume η p > δec −1 . When the average feeding duration of nymphal ticks T = > 0 is small, bi-stability occurs in model (15) in the cooperative feeding and density-dependent grooming case considered in Theorem 1: the trivial equilibrium P 0 and the positive equilibrium (F * + ( ), E * + ( )) are both locally asymptotically stable, and the positive equilibrium We note that in the case considered in Theorem 3, the trivial equilibrium is asymptotically stable since ρ(0) = 0.

Quadri-stability with grooming reactive to tick biting
We consider the case that k = O( −1 ). Consider the expansion k as follows Then the first equation of (10) becomes δη −1 = e −(k 0 +k 1 +o( ))ρ ρ.
For = 0, Equation (31) has two solutions ρ −0 and ρ +0 with 0 < ρ −0 < k −1 0 < ρ +0 if (ek 0 ) −1 > δη −1 . When ρ +0 < r j/2, by solving the quadratic equation (18), we obtain four positive solutions F 0 We now look for For the four positive equilibria and (F −+ ( ), E −+ ( )) when > 0, we suppose that the series expansions are expressed as Using the attachment rate function (18), we get As we did in last subsection, we now use the method of series expansions to look for the signs of the real parts of the corresponding characteristic equation. In what follows, we drop all the subscripts ("−" and "+") for simplicity.
In what follows, we will explore the solutions of H (λ, ) = 0 at each equilibrium The second equation of (10) gives (34) Combining with (32), we have According to the first equation of (10), we can obtain which leads to Substituting (33) into (39), we have Then ρ 1 in (33) can be rewritten as From (36), it is easy to get Based on the obtained expressions of a i j in (22), we get a 11 (λ, ) = −δ + ηe −(k0+k1 +o( ))(ρ0+ρ1 +o( )) (ρ 0 + ρ 1 + o( ))e −(λ0+λ1 +o( ))(τ + ) Substituting a 11 (λ, ), a 21 (λ, ), a 12 (λ, ) and a 22 (λ, ) into the characteristic equation H (λ, ) = 0, we have where we have It is clear that (41) is the characteristic equation for the following systeṁ which is a delay differential equation with positive feedback if δ > A. The stability of (42) is equivalent to the stability of the zero solution for the following ordinary differential equation (see Smith 1987) Then we have the following result: Theorem 5 Assume that δ > A. The stability of equilibria are determined by the signs of A, respectively, i.e., the equilibrium for model (17) is stable (unstable) if and only if A > 0 (A < 0). Now we assume that ηE 0ρ0 e −k 0 ρ 0 < 1.
Proof Using the assumption (45), we can see that the inequality (43) and δ > A are always true. Note that the sign of A is determined by the signs of the following three terms: sign(ρ 0 ), sign(1 − k 0 ρ 0 ) and sign(ηE 0ρ0 e −k 0 ρ 0 ).
Then the sign of A is We remark that in the case considered in the above theorem, ρ(0) = 0 so the trivial equilibrium is always asymptotically stable.

Simulations and discussions
As shown in previous studies (Brauer 1995;Webb 1985;Metz and Diekmann 2014;Magal and Ruan 2018;Wu 1996;Kosovalić et al. 2017Kosovalić et al. , 2014Kosovalić et al. , 2013, algebraicdelay differential systems arise naturally from the population dynamics, when the individuals are physiologically structured and when the total population in a particular physiological stage is subject to certain dynamics leading to an integral equation with the input flow into the stage acting as a continuous forcing.
Here, we formulated a novel coupled system of a differential equation and an algebraic/integral equation to characterize the tick population dynamics, with a particular focus on the tick attachment/fixation behaviors and host grooming response and the implication of this tick-host interaction on the tick population dynamics at the population level. The model formulation is particularly useful to provide insights of multi-stability in different combinations of tick-host pairs. We did consider two special cases: the case of cooperative feeding and density-dependent grooming when bi-stability occurs, and the case of cooperative feeding and grooming in response to tick biting when quadri-stability becomes a feasible scenario.
We now provide a few numerical examples to show bi-stability and quadri-stability are both possible with parameters suggested from experimental settings or field observations (see Table 1).
Bi-stability We start with the case of cooperative feeding and density-dependent grooming, with some parameter values given in the second column of Table 1 and the other parameter values assumed below: A simulation is presented in Fig. 2, from which we observe the existence of three equilibria of model (15): one is the tick-free equilibrium P 0 (0, 0) and the others are positive, the larger positive equilibrium (E * + , F * + ) = (37.1687, 164.9135) is stable and the smaller one (E * − , F * − ) = (5.1798, 7.0223) is unstable. Illustrated also in Fig. 3 are solutions of feeding nymphal ticks with different initial data that approach two stable equilibria P 0 and (E * + , F * + ). Varying the feeding duration T ∈ [0, 9.3718], the two equilibria (in blue) remain to be stable and the middle equilibrium (in red) is unstable, so we have the expected bi-stability and the tick population long-term behaviors depend on the initial condition. If T * = 9.3718, model (15) has only one Fig. 2 Equilibrium of feeding nymphal ticks with respect to average feeding duration T ∈ (6, 9.3718). For each T , model (15) has a stable tick-free equilibrium (bottom blue curve) and two positive equilibria, one stable (top blue curve) and one unstable (middle red curve) unique positive equilibrium (E * , F * ) = (25.9754, 135.9141), and increasing T , the coupled system has the only tick-free equilibrium.
Quadri-stability We provide another set of simulations for the case of cooperative feeding and grooming in response to tick biting. Some parameters are listed in the third column of Table 1 and the remained parameters are fixed as follows: δ = 0.2, k 0 = 2/7, r = 0.02, j = 120.
Due to the incorporation of two time lags (T and τ ) and the integral term in the algebraic equation for F(t), local stability analysis of the multiple equilibria coexisted in the coupled system is very difficult and we are only able to conduct the local stability analysis using the perturbation analysis when the feeding duration T is relatively small comparing with the long life cycle (T + τ ). In the case of cooperative feeding and host grooming in response to tick biting, increasing the feeding duration may lead to nonlinear oscillations around these co-existing equilibria through Hopf bifurcations. In this scenario, obtaining an analytic expression of the critical value of the feeding  (15) with different initial data: (1) when the initial data are fixed at (80, 180), (60, 90), (30, 20) respectively, the numbers of engorged and feeding nymphal ticks approach to (E * + , F * + ) = (37.1687, 164.9135); (2) when the initial data is (5, 8), nymphal ticks will go extinct eventually. Model (15) has two locally asymptotically stable equilibria: (E * + , F * + ) and tick-free equilibrium (0, 0) duration when Hopf bifurcations take place, and describing the global continuation of periodic solutions bifurcated near the corresponding equilibria, remain a challenging task. The global dynamics of such a model has yet to be obtained. Tick population dynamics and tick-borne disease transmission dynamics have been modeled intensively (see, for example, Gaff and Gross 2007;Rosà et al. 2003;Rosà and Pugliese 2007;Wu and Zhang 2020 and references therein). Density-dependent development rates have been empirically estimated using laboratory or field observation, see for example, Ogden et al. (2005). These estimations were used for examining the tick range expansion (Wu et al. 2013). Our contribution here is to use the coupled system to separate the tick attachment/fixation behaviours from the host grooming behaviors in estimating the density-dependent development rates and explore the implication of different combinations of tick attaching and host grooming behaviors. This complements the study of  on tick seeking assumptions and their implications for disease transmission dynamics.
An important step forward is to expand our coupled system for the tick population dynamics to a coupled system for tick-borne disease transmission dynamics when tick population is further stratified by physiological stages and epidemic status. Separating the tick-attaching and host grooming behaviors for ticks at different stages in describing the tick-borne disease transmission dynamics is also important to understand the cooccurrence of ticks at different stages in the same host, which is critically important to understand the role of co-feeding transmission dynamics (Alekseev and Chunikhin  (17) with different initial data: (1) when the initial feeding nymphal ticks are set at 1600, 700, 300, the number of feeding nymphal ticks approaches F −+ = 1443.8375; (2) when the initial feeding nymphal ticks are set at 200, 150, 10, the number of feeding nymphal ticks approaches F +− = 65.1540; (3) when the initial feeding nymphal ticks is 3.3, feeding nymphal ticks will go extinct eventually. Thus, model (17) has three locally asymptotically stable equilibria: (E −+ , F −+ ), (E +− , F +− ) and tick-free equilibrium (0, 0) 1990; Hua et al. 2003;Labuda et al. 1993;Mogl et al. 2011;Ogden et al. 1997;Randolph et al. 2002Randolph et al. , 1996Randolph 2011;Wu and Zhang 2021;Zhang et al. 2017). and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.