Synchronized Tick Population Oscillations Driven by Host Mobility and Spatially Heterogeneous Developmental Delays Combined

We consider a coupled system of delay differential equations for a single-species tick population dynamics, assuming feeding adult ticks are distributed by their hosts in a spatially heterogeneous environment consisting of two patches where egg ticks produced will complete their life cycles with different, normal and diapause, developmental delays. We show that the mobility of adult tick host and the diapause developmental delay combined drive a synchronized oscillation in the total tick populations around a uniquely defined positive equilibrium, and this synchronization makes the oscillatory patterns much simpler in comparison with multi-peak oscillations exhibited in the absence of host mobility.


Introduction
Tick-borne diseases, including tick-borne encephalitis, Lyme disease, anaplasmosis and Crimean-Congo hemorrhagic fever, have emerged as a global public concern since the beginning of the century (Dantas-Torres et al. 2012;Dantas-Torres 2015). In parallel to progress made in the fields of diagnostic and therapeutic technologies (Bologheanu et al. 2020;Serretiello et al. 2020) . 2017), mathematical models and analyses have also been advanced, using differential equations to explore the long-term behaviors of tick-borne pathogen spread in the vector-host interaction (see Wu et al. 2013;Rosà and Pugliese 2007;Zhang and Wu 2020 and references therein).
Tick-borne diseases are mainly transmitted through the bite of ticks. Ticks, second only to mosquitoes for vectoring human, have two families: Ixodidae (hard ticks) and Argasidae (soft ticks) (Hoskins 1991). Approximately 650 Ixodidae species have been reported worldwide (Boulanger et al. 2019). Ticks in the family of Ixodidae distribute widely in natural ecosystems and take major responsibilities for transmitting a variety of tick-borne diseases. Most ticks go through four life stages: egg, larva, nymph and adult. The development in each post-egg stage involves a process of questing, attaching, feeding and engorging (Marquardt 2005), with female adult ticks laying thousands of eggs and dying shortly to complete the life cycle. Due to global warming, tick populations expand their spatial habitats and temporal activity season into winter months (Sagurova et al. 2019). For a better estimation of the tick-borne disease spread, understanding the tick population dynamics is in priority.
An important mechanism for potentially complex tick population dynamics is the developmental heterogeneity due to diapause, an important physiological process for ticks to respond to variable environmental conditions, such as changes in temperature, photoperiod, rainfall and host availability. Tick diapause was described for the first time by Beinarowitch (1907). There are two basic types of tick diapause: behavioral and developmental diapause, which delay host-questing following the transformation of molt and suspend development from one stage to the next, respectively (Dobson et al. 2001;Ogden et al. 2004). Korotkov (2015) estimated the duration of tick development cycle that can reach 3 ∼ 6 years. Shu et al. (2020) proposed a delay differential equation for tick population dynamics with multiple delays and discussed global Hopf branches. Hoch et al. (2010) assumed that all engorged ticks experience diapause and develop until the second spring of next year and investigated the dynamics of tick population with temperature changes. Despite these efforts, quantifying impact of diapause on tick population dynamics remains a challenge partially because the diapause is linked to the experience of ticks in different habitats and environments during their entire life cycle.
Natural habitats have been increasingly separated into many smaller patches by human activities, such as road construction, mining, deforestation and aggressive agricultural cultivation. Although it is unlikely for ticks to move over a long distance, hosts carrying feeding ticks can freely move among different patches. Since the blood meal of hard ticks such as Ixodes ricinus lasts several days (Gray 1981), ticks can be carried over by their hosts and move from one patch to another during the period of blood feeding. Hence, it is highly possible that the patch where a tick quested and attached to a host is different from the patch where the tick drops after blood feeding and engorgement, and hence, the development delays are highly relevant to the patches and their ambient conditions. This calls for patch model formulation of the tick population dynamics.
Some earlier studies have proposed multi-patch population dynamics models, see Arino et al. (2019) and references therein. Guiver et al. (2017) gave an explanation of the effects of dispersal-induced coupling on population persistence across discrete patches. Wei and Wang (2020) described a single-species population model between two patches (the natural reserved and naturally environmental) and revealed the sensitivity of fluctuation intensity and migration rate on the population persistence and extinction. Smith et al. (2014) analyzed the long-term behavior of a structured twopatch metapopulation model and used the perturbation theory to discover an Allee effect for a two-patch model .
Here, taking into consideration of diapause development during some physiological stages of ticks in some specific patches the ticks stay in their life cycle, we propose a tick population dynamics patch model with delay. We incorporate both normal developmental delay and diapause delay in two different environments/patches between which ticks move because of the mobility of respective hosts. We show how tick populations are redistributed between questing and engorgement phases by the mobility patterns of the hosts, and we model the heterogeneity of tick development delays due to their distribution between these two patches, those with normal development and those with diapause delay. For such a coupled system of delay differential equations with multiple delays, we establish the existence and uniqueness of a positive equilibrium and obtain conditions for the local stability. We then conduct some numerical simulations to see how the host mobility and diapause development delay combined synchronize the coupled system to a synchronized mode of oscillation, and how this simplifies the pattern of oscillations in comparison with the situations when the two patches are not connected by host mobility.

Spatial Dynamics of Mobile Hosts
Consider a tick species for which the hosts of feeding ticks in different stages move between two geological patches, dispersing the engorged ticks among these patches characterized by differential development delays between consecutive stages.
For a given host population under consideration, we denote the number of respective host individuals in each of these two locations by H R (t) and H D (t) at time t, respectively, where we use sup-index R for patch facilitating regular development delay and sup-index D for patch facilitating diapause development delay of the ticks in a given developmental stage (larval, nymphal and adult). Ignoring all other population dynamics, the spatial dynamics of the hosts is given by with m h1 and m h2 for the mobility rates of the hosts from patch R to patch D, and from patch D to patch R, respectively. It follows that Based on the method for calculating basic solution matrix of linear ordinary differential equations (Sina and Ali 2003), we have which describes the redistribution dynamics of host individuals (and hence the respective engorged ticks) as follows:

Tick Population Dynamics with Feeding Tick's Mobility Driven by Hosts
We now consider a tick species in a region with two distinct patches where the environmental conditions are suitable for a regular development delay (normally one year) in patch R and a diapause delay (normally two years) in patch D for ticks in each development stage from larva to nymph, and from nymph to adult. We ignore the mobility of ticks during their engorged and questing activity stages, but ticks feeding on hosts move along the hosts which we assume move between the two patches with appropriate mobility rates m l , m n and m a that will further be stratified in terms of the movement from patch R to D or from patch D to R. The flowchart of the evolution of tick population between the two patches with mobility driven by the host migration is given in Fig. 1. Let τ R el and τ D el (τ R el < τ D el ) be the durations which eggs hatch into questing larval ticks and then become feeding larval ticks for the first blood meal on the patch R and patch D, respectively. ρ R e f l and ρ D e f l are corresponding probabilities of survivability from eggs to feeding larval ticks in the patch R and patch D, respectively. The dynamics of feeding larval ticks takes the following form . (1) Since blood feeding usually takes several days, larval ticks are able to be carried over by their hosts during this period. Let τ f l be the feeding time of larval ticks and m l be Fig. 1 A schematic illustration of tick population dynamics with mobility driven by host migration between patch R and patch D: female producing adult ticks after breeding on a host drop to the ground to lay eggs, and then, eggs hatch into the larval ticks. Larvae (L q ) start to quest their first host for an initial blood meal. After feeding, larval ticks (L f ) become engorged larval ticks (L e ) and molt into nymphs. Nymphs experience the same procedure, questing (N q ), feeding (N f ) and engorged (N e ), for a second blood meal, and eventually molt into adult ticks. The adult ticks ( A q ) then hunt a larger host, such as deer, where they are able to feed (A f ) and become engorged (A e ), and finally drop from host for the reproduction ( A p ). There are three mobility patterns, characterized by the random movements of host for feeding larvae, nymphs and adult ticks with migration rates m l , m n and m a , respectively the migration probability of host population between the two patches, then we have After the engorged larval ticks fall to the ground, they develop into nymphs with the survival probabilities ρ R e f n and ρ D e f n , respectively. The corresponding developmental delays are τ R en and τ D en (τ R en < τ D en ) on each patch. The evolution can be described as follows .
Similarly, let m n and m a be the migration probabilities of host population which can carry feeding nymphal ticks and feeding adult ticks and redistribute them in patch R and patch D when feeding ticks become engorged and drop off from the hosts and τ f n and τ f a be the feeding durations of nymphal ticks and adult ticks, respectively. Let ρ R e f a and ρ D e f a be survival probabilities of adult ticks developing from engorged nymphs and τ R ea and τ D ea be the corresponding development time delays. Then, we yield , Following engorgement, adult ticks leave from host and start to be ready for their reproduction. Here, let ρ R epa and ρ D epa be development rates from engorged adults to egg-producing adults, and the corresponding development time delays are τ R epa and τ D epa , respectively. We describe this process by .
( 2) Finally, these egg-producing adult ticks will reproduce by laying eggs and die shortly after that. The dynamics of eggs can be written as where f R and f D are reproduction functions on patch R and patch D and both use with Ricker functions, that is, where p R and p D are maximal numbers of eggs that an egg-producing female adult tick can lay per unit time, s R and s D represent the strengths of density dependence, and μ R and μ D are exit rates of eggs including mortality rate and developmental rate from egg to larval stage. From Eqs.

The Case of Distinct Mobility of Hosts for Feeding Adult Ticks
An important case is when regular development takes place in patch R, while diapause development occurs in patch D. Let the corresponding time delays be τ R and τ D for the entire life cycle, respectively. This naturally leads to our standing assumption in this study: τ D = 2τ R . Using τ = τ R , we have τ D = 2τ . Since the survival probability in patch D is smaller comparing with that in patch R, we write with κ l , κ n , κ a , κ p ∈ (0, 1). We assume two patches do not have any distinction for the birth characteristics (the maximum reproduction numbers of eggs, the exit rates from the egg stage and the strengths of density dependence), leading to We do not focus on the case of insignificant mobilities of hosts for larval and nymphal ticks, but assume there is a preference of hosts for feeding adult ticks moving from patch R to patch D. Therefore, we have wherem a is the average mobility rate of the hosts for feeding adult ticks between the two patches and 0 < δ < 1. We can now drop all the subscripts of the parameters in (5) and (6), and we can also calculate all entries in matrix Q as follows: where ρ = ρ R epa ρ R e f a ρ R e f n ρ R e f l , Δ = 2m a τ f a .
Similarly, for J (y 2 ), we get

Fig. 3 A schematic illustration of the function R(y 1 )
Therefore, Comparing with Thus, there exists a unique positive equilibrium (y * 1 , y * 2 ) for (10). Based on the analysis for J (y 2 ) and R(y 1 ), we have the following conclusions about the location of the unique positive equilibrium (y * 1 , y * 2 ): (iii) if min{ p }, without loss of generality, we may assume that u 22 < u 11 . Figure 4c describes the schematic diagram of J (y 2 ) and R(y 1 ).
In summary, we have the following result: Theorem 1 There exists a unique positive equilibrium (E R * , E D * ) for model (4).

Stability
The linearized system (4) at the tick-free equilibrium P 0 = (0, 0) is given by This gives a linear system of delay differential equations with positive feedback as the coefficient matrix p q 11 q 12 q 21 q 22 > 0.
Therefore, the stability of the zero solution of (12) is equivalent to the stability of the zero solution for the corresponding linear ordinary differential equation The characteristic equation of model (13) at the zero equilibrium is Moreover, we have Then, the two eigenvalues of (13) are Clearly, the two eigenvalues are both negative if the following condition is satisfied Thus, we have established the following result: Theorem 2 Assume (14) holds. Then the trivial equilibrium (0, 0) of model (4) is locally asymptotically stable.
In what follows, we consider the stability of the positive equilibrium P * of model (4). Translating the positive equilibrium into the origin through introducing (4) and linearizing it, we have where We can further compute It is clear that (15) is a positive feedback system if the following condition is satisfied From (8) and (11), the condition (16) can be guaranteed if which are equivalent to 2k p k a k n k l ρ .
Then, the stability of the positive equilibrium (E R * , E D * ) can be determined by the stability of the zero solution to the following linear system of ODE equations: The characteristic equation of (18) is Its roots are Since (17) implies γ 1 and γ 2 are both positive, the stability of the positive equilibrium (19) can be rewritten as μ 2 + γ 1 γ 2 (q 11 q 22 − q 12 q 21 ) − μ(γ 1 q 11 + γ 2 q 22 ) > 0, which is equivalent to It is easy to see that (20) is satisfied if the following condition holds Therefore, we have Theorem 3 Assume that (17) and (21) hold. Then, the positive equilibrium P * (E R * , E D * ) is locally stable.
We now consider the case when the condition of a positive feedback system is not satisfied. The characterization equation of (15) is From the analysis above, we can see that when τ = 0, the positive equilibrium P * is locally asymptotically stable if (21) is satisfied. In the absence of a positive feedback at the non-trivial equilibrium, Hopf bifurcation can take place when the characteristic equation has a pair of purely imaginary zeros. We now discuss whether (22) has a pair of purely imaginary roots.
Since γ 1 and γ 2 are both functions of E R * and E D * , it is difficult to solve the transcendental equation (24) in closed forms. However, given our proof of the unique positive equilibrium, all relevant coefficients are specifically given once model parameters are given. Therefore, we can develop a continuation procedure to locate the minimal value of τ , for a given mobility δ ∈ [0, 1], where (24) has a solution ω > 0. In particular, sincem a = 0 corresponds to the situation when two patches are isolated from each other, there exists a critical value τ R * of delay τ such that the population starts to oscillate in patch R if τ > τ R * , and there exists a critical value τ D * of the delay such that the population starts to oscillate in patch D if τ > τ D * in patch D and if τ > τ * D (see Zhang and Wu (2019) for these critical values in closed forms). Normally, τ R * < τ D * . We will concentrate the case where there is a positive mobilitym a and the delay is within [τ R * , τ D * ], so tick population in patch R alone will converge to the positive equilibrium, while tick populations in patch D exhibit oscillatory patterns as a local Hopf bifurcation of periodic solutions takes place near the positive equilibrium.
Whenm a starts to increase from 0, two patches are connected by host mobility and system (31) becomes a transcendental function involving sin(nωτ ) and cos(nωτ ) with n = 1, 2, 3. We look for solutions of (ω, τ ) as a continuous function from (ω D * , τ D * ) with δ as the varying parameter. However, a combination of the lack of a close form of the equilibrium as a function of the parameters including the delay and the mobility rate, the model as a system rather than a scalar equation and the existence of multiple delays makes the Hopf bifurcation analysis so complicated that the analytic calculation for the critical value and the Lyapunov coefficient of a Hopf bifurcation of periodic solution is not feasible. We will design some numerical simulations in next section based on the estimated critical values of delay and mobility from a continuation procedure.

Numerical Simulation
In this study, we consider the homogeneous situation, so seasonal variations are incorporated by taking the unit time as one year to normalize the parameter values taken from some published literature. We refer to the monograph Zhang and Wu (2020) for the original sources of the parameter values: We also fix all proportions of development rates from one stage to the next at the same value κ, i.e., κ = κ p = κ a = κ n = κ l = 0.55.
Without mobility, the model (4) has a positive equilibrium E * (14970, 65090). Considering the delay τ as the bifurcation parameter, we can calculate the critical value τ * = 0.69. From Fig. 5a, we observe that the egg density in regular patch E R (t) is locally asymptotically stable when τ = 0.55. The qualitative behavior of egg density E R (t) in model (4) undergoes periodic solution, dual-peak and more complex oscillation with increasing delay τ , which is illustrated in Fig. 5: (b) τ = 0.7; (c) τ = 4; (d) τ = 10. However, the solutions of egg density in the diapause patch E D (t) do not change with different delays and remain stable. Note that when we increase the proportion κ into 0.7, Fig. 6 compared with Fig. 5c and d shows that oscillation of egg density in the diapause patch starts to emerge. This shows that lower survival probability can keep egg dynamics in diapause patch to be stable.
When the host movement is considered between patch R and patch D, and using the mobilitym a = 146 as an illustration, we show that feeding adult ticks will be carried over between the two patches and the positive equilibrium becomes (37,550,11,620). Therefore, first of all, the host mobility sustains the persistence of tick population in the both patches, though the regular patch is the preferred habitat. This is illustrated in Fig. 7a and b.
With the same value of delays as in Fig. 5, we describe the dynamics of model (4). From Fig. 7b, it is clear that the egg density is stable when τ = 0.7 and the critical value of time delay of model (4) with host movement is getting larger, τ * * = 2.12. With increasing value of the time delay, the behaviors of egg densities in both patches  Fig. 5c and d, solutions of egg density for the model (4) for κ = 0.7. Clearly, increasing κ will contribute the oscillation phenomenon for egg density in the diapause patch  (4) with host mobility between the two patches R and D for τ = 4, where blue line denotes the nymph density in regular patch and the red line represents the one in diapause patch; other parameters are the same as those in Fig. 7 (Color figure online) become periodic, as shown in Fig. 7c and d. We also plot, when τ = 4, the engorged nymph densities in Fig. 8. We also compare Fig. 7b-d with Fig. 5b-d, and we note that egg densities in the regular patch are asymptotically constant with host mobilities when τ = 0.7 and exhibit much simpler single-peak oscillatory pattern when τ = 4 and τ = 10 in comparison with dual-peak oscillations and complex oscillation in Fig.  5b-d. Hence, host mobility not only synchronizes but also simplifies the collective behaviors.
In order to better understand the parameter values listed above, we normalize the ) and obtain whereτ R = 365τ R andτ D = 365τ D ;Ẽ R andẼ D represent the total numbers of eggs to be hatched within one day. Based on the process of normalization, it means that tick population will take 365τ R or 365τ D days in each patch to complete their life cycle from the view of biology. For instance, τ = 0.55 in Fig. 5a is equivalent toτ = 0.55 × 365 = 200.75, which implies that the life cycle of tick population in regular patch is 200.75 days.

Discussion
It is commonly known that delay, if increased passing a critical value, in a negative feedback system destabilizes the system and generates nonlinear oscillations through the Hopf bifurcation mechanism. It is also known that increasing the delay further may lead to complicated oscillatory patterns including multi-peaks with a single period. Multiple and large delays may appear in tick population dynamics due to the socalled diapause development, when ticks suspend their development in response to unfavorable conditions. It was established in Shu et al. (2020) and Zhang and Wu (2019) that in single spatial location, when a portion of ticks undergo diapause during their life cycle, nonlinear oscillations can take place.
Here, we consider a landscape consisting of two patches (patch R and patch D) where ticks in different patches will experience regular (in patch R) and diapause (in patch D) development. When these are in isolation in terms of population dynamics of a particular tick species, we may observe patterns of stabilization to a positive equilibrium in one patch and oscillation in another patch with potentially multi-peaks within a given period. We then show the mobility of hosts on which the adult ticks feed for development can redistribute engorged adult ticks in two patches, creating two cohorts of ticks developing in two different habitats until they feed on hosts during the questing activities in the adult stage. The tick population dynamics is then described by a couple system of delay differential equations with two delays, the normal and diapause development delays. We show that the model system may have a unique positive equilibrium and its stability may change due to the mobility of the hosts and the diapause delay in combination. This will tend to synchronize the system to periodic oscillations with single peak within a given period. This synchronization to a simple-looking periodic solution seems to be new and may offer another angle to view tick population dynamics in the natural environment. It remains for further studies to see how seasonal variation of the habitat conditions coupled with the diapause delay can impact the tick population dynamics (this requires a periodic system of delay differential equations with time-dependent delays) and influence the patterns of tick-borne disease spread.