Discrete stage-structured tick population dynamical system with diapause and control

1.   Introduction

Diapause is a key survival mechanism of Ixodes and other invertebrates, such as mosquitoes, dragonflies and ladybugs ^[1,2,3], which synchronizes the rhythm of the life cycle with that of the environment to ensure that the active phase of the life cycle occurs when food is abundant and other aspects of the environment are conducive to survival ^[4]. Diapause is also a developmental stagnation period caused by adverse climatic conditions, particularly photoperiod and relative humidity ^[5].

In the natural world, many species go through some of the distinct life stages, and while individuals at each stage are biologically identical, subgroups of these species differ in physical characteristics and have different vital behaviors. The single-species population dynamic model with stage-structured^[6,7,8] has attracted the attention of many scholars. Ticks, as vectors of Lyme disease, tick-borne encephalitis and human babesiosis, are vectors of vector-borne diseases that have a significant impact on human health. Since ticks respond differently to the environment at various stage, The life cycle of ticks consists of four successive developmental stages, namely egg, larva, nymph and adult ^[9,10]. The ticks has different hosts at different stages. For larvae and nymphs, rodents are the most important reservoirs, such as the white-footed mouse. The most important hosts of adults is larger mammals such as ungulates, the most common white-tailed deers ^[11]. After the eggs evolve into larvae, and the larvae develop into nymphs through a blood meal on a rodent host. The nymphs continue to seek the host and feed. After engorging, the nymphs will drop off the host and go through a period of development, after which they evolve into adults with the male continuing to seek a host while the female lays eggs. The diapause of the tick population is diverse and can occur at any stage of their life cycle. The basic diapause types can be divided into two categories: developmental diapause (temporary suspension of engorged tick development) and behavioral diapause (interruption of host-seeking activity of unfed ticks). Belozerov et al. ^[12] studied the existence of photoperiod controlled diapause in the nymphs of prostriate ticks and its influence on the nymphs dynamics through some data and confirmed that behavioral diapause in adults of other tick species was a well established phenomenon ^[13]. Gray ^[14] found that the diapause behavior was terminated when the female tick approached the diapause male tick, demonstrating a richer dynamic behavior of ticks.

In recent years, many scholars have considered diapause in the life cycle of population and studied the important influence of diapause on population growth and development by establishing some reasonable mathematical models. Lou et al. ^[15] regarded the diapause period as a dynamic process independent of the normal period, and studied the influence of diapause on mosquito population dynamics by constructing differential equation model. The sensitivity analysis of the parameters related to diapause proved that by reducing the diapause mortality, the short diapause period could increase the survival ability of mosquitoes, which was more conducive to the normal growth of subsequent mosquito populations. Zhang and Wu ^[16] analyzed a population dynamics model with ticks as vectors by combining the development delay and diapause delay and computing the Hopf bifurcation value by introducing and analyzing the parametric trigonometric functions. Shu ^[17] analyzed a time-delay differential equation for a diapause tick population and performed a Hopf bifurcation analysis of the model, showing rich tick dynamic behavior. However, mathematical models considering the role of diapause in the life cycle of ticks is still rare, and most of the existing studies are based on differential equation dynamical systems, ignoring the critical feature that there is no generational overlap in tick populations. Consequently, it is more practical to consider a discrete dynamic model.

At present, the main measures for tick control include spraying acaricides, using tick-repellent and vaccination. Existing vaccines are still in the stage of screening and validation of tick functional molecules to find suitable antigens. Although there is a small number of commercialization vaccines were developed, such as Boophilus microplus Bm86 and Bm95 antigen vaccines (TickGARD/Gavac), but were not widely available due to their inconsistent effectiveness ^[18]. Acaricide spraying has been a relatively effective way to control tick populations in recent years. Quantifying and incorporating control measures into tick populations dynamics models will help us more comprehensive explore the development process of ticks and provide theoretical reference for controlling tick-borne disease transmission.

Our paper is organized as follows. In Section 2, we combine the impact of diapause on the life cycle of tick populations, divide immature ticks into diapause immature and normal immature, and construct a stage-structured non-monotonic tick population dynamical system. In Section 3, we calculate the net reproduction number and discuss the global stability of tick-free fixed point; we also confirm the existence of a unique coexistence fixed point, and infer its local stability under some specific conditions. Moreover, we conclude that when the net reproduction number is at the critical value, the model will undergo a transcritical bifurcation. In Sections 4 and 5, we discuss the uniform persistence and investigate the effects of different acaricides spraying strategies on tick population dynamics, respectively. In Section 6, some numerical simulations are given to support the theoretical analysis and exhibit these rich tick population dynamical behaviors. We end the paper with a brief discussion.

2.   The model formulation

Denote $I(t), I_d(t),$ and $M(t)$ as the amounts of normal immature tick population, diapause immature tick population and mature tick population at time t, respectively. We have the following diagram depicting the growth of tick populations(see Figure 1). Let $b$ be the oviposition rate; $\gamma$ denotes the survival rate of eggs; $s_0$ is the hatching probability from eggs to the larval stage; $d$ is the diapause rate of immature ticks; $\beta$ represents the transition rate from immature to mature(i.e., transition rate from nymphs to adults); $\delta$ represents the exit rate of diapause immature ticks. We assume that normal and diapause tick population in immature competes for each other and independent of the mature stage. $S_1(x)$ is the Beverton-Holt type $\frac{c}{a+x}$ ^[19,20] nonlinear function to describe the survival rate of immature, and we assume that survival rates of diapause immature ticks and mature ticks are constant $S_d$ and $S_m,$ respectively. We let $E_m$ represents the acaricidal effect of spraying acaricides. Following the block diagram and assumptions described above, and we consider the following difference system:

Since $S_1(I(t)+I_d(t)) = \frac{c}{a+I(t)+I_d(t)},$ it satisfies the following conditions:

For the view of biology, we assume that $\delta+d > 1.$

Denote $z(t)$ : = $(I(t), I_d(t), M(t))^{T}$. The following matrix form is used to represent model $(2.1):$

where the projection matrix $\Phi(z)$ is given by

Let $A = (a_{ij})\in \mathbb{R}^{n\times m}$ and $B = (b_{ij})\in \mathbb{R}^{n\times m}.$ Define $A\leq B$ if and only if $a_{ij}\leq b_{ij},$ for all $i = 1, ..., n,$ and $j = 1, ..., m.$ Therefore, from $(\Sigma 1)$ we can obtain that $\Phi(x)\geq \Phi(y),$ for any $x\leq y.$

Furthermore, combined with conditions $(\Sigma 1),$ we have

Let $\widetilde{M}(t+1)$ satisfy the recursion $\widetilde{M}(t+1) = \beta c+S_m(1-E_m)\widetilde{M}(t).$ Then we can easily deduce that $\widetilde{M}(t+1) = \beta c\sum\limits_{i = 0}^t(S_m(1-E_m))^{i}+(S_m(1-E_m))^{t+1}\widetilde{M}(0)$ and $\widetilde{M}(t)$ has limit at $t\rightarrow \infty.$ Since $0\leq M(t)\leq \widetilde{M}(t),$ $M(t)$ is bounded as $t\rightarrow \infty.$ Let $\overline{M}$ be its upper bound, that is, $M(t)\leq \overline{M}, \; t = 1, 2, ...$. Similarly, we can prove that $\overline{I_d}$ be the upper bound of $I_d(t)$ and $I_d(t)\leq \overline{I_d}, \; t = 1, 2, ...$. Then, we can follow from model $(2.1)$ that

Therefore, the positivity and boundedness of the solutions of the model $(2.1)$ can be summarized as follows

Theorem 1. System $(2.2)$ is point dissipative. Let $(I, I_d, M)$ be the solution of model $(2.1).$ Denote

Then $\Upsilon$ is positively invariant under the flows of the system $(2.2)$ and is attracting to all solutions of $(2.2)$ under condition $(\Sigma 1).$ That is to say, there is a compact set $\Upsilon \in \mathbb{R}_{+}^3$ such that every forward solution sequence of the system $(2.2)$ enters $\Upsilon$ in at most two-time steps, and remains in $\Upsilon$ forever after.

3.   The net reproduction number

We define the net reproduction number $R_0$ of the tick population using the methods in ^[21,22,23], which is the expected amount of descendants produced by an individual over the course of its life.

We can obtain the following fertility matrix

and transition matrix

Then the inherent projection matrix is $\Phi(0) = F+T(0)$, and the next generation matrix is given by

where $\mathcal{I}$ denotes the identity matrix. Solving $\det(\lambda \mathcal{I}-F(\mathcal{I}-T(0))^{-1}) = 0$, we can yield the net reproduction number

It is clear that the tick-free fixed point $E_0(0, 0, 0)$ always exists for the model $(2.1).$ For any $z(0)\geq 0,$ we have $z(t)\geq 0,$ it follows from system $(2.2)$ and condition $(\Sigma 1)$ that

and

and repeating this manner, we can get that

It is easily verify that the inherent projection matrix $\Phi(0)$ is non-negative, irreducible, and primitive. Thus, $\Phi(0)$ has a simple, positive and strictly dominant eigenvalue $r.$ If $R_0 < 1,$ we can obtain $r < 1$ and

by using the methods in ^[24]. It indicates that $E_0(0, 0, 0)$ is globally asymptotically stable. In addition, when $R_0 > 1,$ $\Phi(0)$ has a positive eigenvalue larger than one. Thus, we linearize model $(2.1)$ at the tick-free fixed point $E_0(0, 0, 0)$ and obtain a positive eigenvalue larger than one. Clearly, the tick-free fixed point $E_0(0, 0, 0)$ is unstable. Therefore, We have the following stability result of the tick-free fixed point $E_0.$

Theorem 2. If $R_0 < 1,$ model $(2.1)$ only has a unique tick-free fixed point $E_0(0, 0, 0),$ it is globally asymptotically stable, and is unstable if $R_0 > 1.$

3.1.   Transcritical bifurcation

We define $f: = \Phi(z(t))z(t)$ be the map from $\mathbb{R}_{+}^{3} \to \mathbb{R}_{+}^{3}.$ Then we linearize the model $(2.1)$ at $E_0(0, 0, 0)$ and obtain the relevant Jacobian matrix

We can obtain that the characteristic equation of $(3.2)$ is

where

Combined with the characteristic equation $(3.3)$, we can find that $f_1(1) = 1+b_1+b_2+b_3 = 0$ is equivalent to $R_0 = 1,$ and $f_1(\lambda) = (\lambda-1)(\lambda^2-(b_2+b_3)\lambda-b_3) = 0.$ Therefore, the other two eigenvalues of Jacobian matrix $(3.2)$

It is easy to derive that when $R_0 = 1$ and $-1 < b_2 < min\{2-b_3, 1-2b_3\}$ or $-1 < b_3 < 0,$ the eigenvalues of Jacobian matrix $(3.2)$ $|\lambda_{1, 2}| < 1, \lambda_3 = 1.$ Therefore, we have the following theorem:

Theorem 3. If $R_0 = 1,$ and $-1 < b_2 < min\{2-b_3, 1-2b_3\}$ or $-1 < b_3 < 0,$ model $(2.1)$ undergoes a transcritical bifurcation.

Proof. We use $b$ as the bifurcation parameter, and then $R_0 = 1$ is equivalent to

Let $B(t) = b-b^*$ as a new dependent variable. Expanding $(2.1)$ as a Taylor series at $(I(t), I_d(t), M(t), B(t)) = (0, 0, 0, 0)$ up to terms of 2 produces

where $\mathcal{F} = \left(\begin{array}{cccccccccc} f_1(I(t), I_d(t), M(t), B(t)) \\ f_2(I(t), I_d(t), M(t), B(t)) \\ f_3(I(t), I_d(t), M(t), B(t)) \\ 0 \end{array}\right)$ is shown in Appendix.

We construct an invertible matrix

Under the transformation

system $(3.4)$ becomes

where $\widetilde{\mathcal{F}} = \left(\begin{array}{cccccccccc} \widetilde{f_1}(u(t), v(t), w(t), \mu(t)) \\ \widetilde{f_2}(u(t), v(t), w(t), \mu(t)) \\ \widetilde{f_3}(u(t), v(t), w(t), \mu(t)) \\ 0 \end{array}\right)$ is shown in Appendix.

By the center manifold theory ^[25], the stability of $(u(t), v(t), w(t)) = (0, 0, 0)$ near $\mu(t) = 0$ can be determined by study a family of parameters equations on a center manifold which can be expressed as

where $w(t)$ and $\mu(t)$ are sufficiently small. In order to compute the center manifold and determine the equation on the center manifold, we assume

where $o(|(w(t), \mu(t))|^2)$ means terms of order greater than 2 in the combination of $(w(t), \mu(t)).$ Substituting $(3.6)$ into the center manifold equation, we have

By simple computation for center manifold, we obtain

Thus, model $(3.5)$ restricted to the center manifold is shown below

It is easy to obtain that $(3.8)$ satisfies the following conditions

Therefore, it can be concluded that $(2.1)$ undergoes a transcritical bifurcation ^[25,26,27].

3.2.   The coexistence fixed point

If the system $(2.2)$ exists a coexistence fixed point $E^{*}(I^{*}, I_d^{*}, M^{*}),$ then its components must fit in with:

From the second and third equation of $(3.9),$ we can conclude that

Substituting $(3.10)$ into the first equation in $(3.9),$ we can obtain

Thus

where

Since $a_1 = \frac{c}{a},$ we can yield

Note that $I+I_d = 0$ implies $M = 0.$ We find that if and only if $I+I_d > 0,$ $a_1\kappa-(1-S_m(1-E_m))(1-(1-\delta)S_d) > 0,$ namely, $R_0 > 1.$ From $(3.10)$ and $(3.12),$ we can compute unique coexistence fixed point $E^*$ of model $(2.1).$

Theorem 4. If and only if $R_0 > 1,$ model $(2.1)$ exists a unique coexistence fixed point $E^{*}(I^{*}, I_d^{*}, M^{*}).$

Assume $R_0 > 1,$ and let $W(t) = I(t)+I_d(t),$ through simple calculation, we have

For simplicity, we define $S_1^{'}(W)$ as $\frac{\partial S_1(W)}{\partial I}$ and $\frac{\partial S_1(W)}{\partial I_d},$ define $(S_1(W)I)^{'}$ as $\frac{\partial S_1(W)I}{\partial I}.$ Then the Jacobian matrix of model $(2.1)$ at $E^{*}$ can be given as follows

The characteristic polynomial of $J_f(E^{*})$ is given by

where

Theorem 5. When $R_0 > 1,$ the coexistence fixed point $E^{*}(I^{*}, I_d^{*}, M^{*})$ is locally asymptotically stable if

Proof. From $(3.11),$ we can conclude that

Substituting $(3.15)$ into $f_{2}(1) = 1+c_1+c_2+c_3,$ and combining with condition $(\Sigma 1),$ we can get

Note that

Then if $c_1+c_3 < 0,$ we can get $-f_{2}(-1) > f_{2}(1) > 0.$ Based on previous assumptions about parameters, we have $c_3 > 0,$ then $c_1 < 0, \; \; \; \lvert c_1 \rvert > \lvert c_3 \rvert,$ and $1+c_2 > f_{2}(1) > 0.$ Thus, we have $c_2-c_1c_3+1-c_3^2 = 1+c_2-c_3(c_1+c_3) > 0,$ namely,

From $c_1 < 0,$ we have $-c_1 = (1-\beta)S_1^{'}(W)I+(1-d)(1-\beta)S_1(W)+(1-\delta)S_d+S_m(1-E_m) > 0,$ which is equivalent to

Thus,

That is to say, $f_{2}(1) < 1$ and $-c_2 > c_1+c_3 > c_1.$ Then we get

Therefore, if $c_3-c_1 < 1,$ then we have

Combining with $(3.16),$ we have $\lvert c_2-c_1c_3 \rvert < 1-c_3^2.$ From Jury criterion ^[28,29], it follows that the roots of the characteristic polynomial $f_2(\lambda)$ have less than one. Thus, $E^{*}(I^{*}, I_d^{*}, M^{*})$ is locally asymptotically stable.

4.   Uniform persistence

It follows from Theorem $2$ and $4$ that when $R_0 > 1,$ the tick-free fixed point $E_0$ is unstable and there is a unique coexistence fixed point. Let $f: = \Phi(z(t))z(t)$ be the map from $\mathbb{R}_+^{3}$ to $\mathbb{R}_+^{3},$ and take $\mathcal{H}$ as the boundary of $\Upsilon$ defined in $(2.3).$ Clearly, we can conclude that $f^{t}(\Upsilon \backslash \mathcal{H}) \subset \Upsilon \backslash \mathcal{H}$ from Theorem $1,$ where $f^{t}(x)$ represents the $t$th iteration of $x$ under the map $f.$ By Theorem 2.1 in ^[24], it indicates that there exists a global attractor $\Lambda$ in $\Upsilon.$ Let $N: = \{(0, 0, 0)\} \subset \mathcal{H}$ be the largest compact invariant set, then $\Upsilon \backslash \; N$ is positively invariant.

Since $\Phi(0)$ is non-negative and irreducible, it has a dominant positive eigenvalue that we denote it as $r > 1,$ and has a relevant positive left eigenvector $q > 0,$ such that

Choose $r^{*}\in(1, r).$ Then

Define a vector norm $\Vert \cdot \Vert$ such that $\Vert z \Vert: = q^{'}z.$ By the continuity of $\Phi(z),$ there exists $\rho > 0,$ such that

for all $z\in U: = \{z\in \mathbb{R}_{+}^{3}, \; \; \Vert z \Vert \leq \rho\}.$ Let $z(t)$ be a solution of system $(2.2).$ If $z(t)\in U,$ it follows from $(4)$ that

which is equivalent to $\Vert z(t+1) \Vert > r^{*} \Vert z(t) \Vert > \Vert z(t) \Vert,$ for all $z\in U.$ Thus, $\lim\limits_{t\to \infty}\inf z(t)\geq \rho$ for all non-zero orbits in $\mathbb{R}_+^{3},$ which indicates that the tick-free fixed point is a uniform repeller. By Theorem 2.1 in ^[30], it means that

$(1)$ $N$ is isolated in $\Upsilon,$

$(2)$ $S(N)\subset N,$ where $S(N)$ represents the set of points whose solution sequence for system $(2.2)$ converges to $N,$ which implies that

$(1)$ $N$ is isolated in $\Lambda,$

$(2)$ $S(N)\subset \mathcal{H}.$

Therefore, system $(2.2)$ is uniformly persistent. The results can be summarized as follows

Theorem 6. If $R_0 > 1,$ system $(2.2)$ is uniformly persistent.

5.   Acaricidal effect of spraying acaricides

In this section, we will mainly concentrate the impact of spraying acaricides on tick population control. For simplicity, we ignore the diapause phenomenon in the model $(2.1).$ Next, we will consider two patterns of spraying acaricides: constant spraying and periodic spraying.

5.1.   Constant spraying

Model $(2.1)$ is simplified to the following form

$S_1(I)$ also satisfies $(\Sigma 1).$ The matrix form of model $(5.1)$ can be given as follows:

Note that $\Phi_{cs}(z_{cs})$ and $\Phi(z)$ have the same properties. It is obvious that model $(5.1)$ exists a tick-free fixed point $E_{cs}(0, 0),$ and similar to the proof of Theorem 2, we can obtain its global asymptotic stability. If the model $(5.1)$ exists a coexistence fixed point $E_{cs}^*(I_{cs}^*, M_{cs}^*),$ it satisfies

which implies that $I$ is a solution of the following equation

Hence, we can conclude that $(5.1)$ exists a coexistence fixed point $E_{cs}^*(I_{cs}^*, M_{cs}^*)$ if and only if

Similarly, we can prove that model $(5.1)$ is also point dissipative. We know that $1\leq M(t)\leq \overline{M}$ from Theorem 1, then $0\leq I(t)\leq bs_0\gamma S_m(1-E_m)\overline{M}+(1-\beta)c.$ Thus, we can obtain a compact invariant set

such that every forward solution sequence of $(5.1)$ enters $\Upsilon_{cs}$ in at most two-time steps, and remains in $\Upsilon_{cs}$ forever after.

Theorem 7. If $R_{cs}^0 > 1,$ then the unique coexistence fixed point $E_{cs}^*(I_{cs}^*, M_{cs}^*)$ of model $(5.1)$ is globally asymptotically stable.

Proof. Define the map $f_{cs}: \mathbb{R}_{+}^{2} \to \mathbb{R}_{+}^{2}$ for the right-hand of model $(5.1).$ The Jacobian matrix of model $(5.1)$ at $E_{cs}^*(I_{cs}^*, M_{cs}^*)$ yields the following form

The characteristic equation of $J_{f_{cs}}(E_{cs}^*)$ is

where

$\text{tr}(J_{f_{cs}}(E_{cs}^*)) = (1-\beta)(S_1(I_{cs}^*)I_{cs}^*)'+S_m(1-E_m),$

$\text{det}(J_{f_{cs}}(E_{cs}^*)) = (1-\beta)(S_1(I_{cs}^*)I_{cs}^*)'S_m(1-E_m)-bs_0\gamma S_m(1-E_m)\beta(S_1(I_{cs}^*)I_{cs}^*)'.$

Combining with Eq $(5.2),$ we can conclude that

and

Therefore, it follows from ^[31] that $E_{cs}^*(I_{cs}^*, M_{cs}^*)$ is locally asymptotically stable.

Since $J_{f_{cs}}(z_{cs})$ is a non-negative matrix for all $z_{cs},$ then $f_{cs}(z_{cs})$ is monotone. It is worth noticing that every solution starting on the boundary of $\mathbb{R}_{+}^{2}$ except $(0, 0)$ enters the positively invariant set int($\mathbb{R}_{+}^{2}$) in at most two-time steps. Now, we choose $z_{cs}(0) = (I(0), M(0))\in \mathbb{R}_{+}^{2}.$ According to the definition of the compact set $\Upsilon_{cs},$ then we can conclude $z_{cs}(0)\in \text{int}(\mathbb{R}_+^{2})\cap\Upsilon_{cs}.$ Obviously, $E_{cs}^*\in \Upsilon_{cs}.$ Let $\bar{\xi}$ be the largest element in $\Upsilon_{cs},$ namely, $\bar\xi = \sup\Upsilon_{cs} = (bs_0 \gamma S_m(1-E_m)\overline{M}+(1-\beta)c, \overline{M}),$ then we can obtain $f_{cs}(\bar{\xi})\leq \bar{\xi}.$ Due to $\Phi_{cs}(0)$ is an irreducible non-negative matrix, the spectral radius $r(>1)$ of $\Phi_{cs}(0)$ with its relevant positive eigenvector $\upsilon$ such that $\Phi_{cs}(0)\upsilon = r\upsilon.$ If $r > 1,$ we can obtain that $f_{cs}(\epsilon\upsilon) = r\epsilon\upsilon+o(\epsilon)\geq\epsilon\upsilon$($\epsilon$ is small enough). Therefore, for given $z_{cs}(0)\in \text{int}(\mathbb{R}_+^{2})\cap\Upsilon_{cs},$ by defining $\epsilon$ small enough that $\underline \xi \equiv \epsilon\upsilon\leq z_{cs}(0)$ and $\underline \xi \leq f_{cs}(\underline \xi).$ Therefore, we can conclude from ^[28] that $E_{cs}^{*}$ is globally attractive. Clearly, $E_{cs}^{*}$ is globally asymptotically stable.

5.2.   Periodic spraying

We assume $E_m = E_m(t)$ is a 2-periodic function in the model $(5.1),$ then model $(5.1)$ becomes

We consider half a year as a unit of time, and assume $E_m(0) = E_1 > 0, E_m(1) = E_2 > 0, E_m(2) = E_1, E_m(3) = E_2, ....$ The methods in ^[22] can be used to construct the critical value $R_{ps}^0$. A periodic linear system can be considered as follows

where

Through the definition of $E_m(t),$ we can obatin

where

Therefore, the projection matrix over a full cycle consisting of two time units is given by the following matrix

where

Similarly, we can compute $R_{ps}^0,$ which is the positive strictly dominant eigenvalue of the matrix $F_{ps}(\mathcal{I}-T_{ps})^{-1},$ and its form is as follows

where

The characteristic equation of $F_{ps}(\mathcal{I}-T_{ps})^{-1}$ is

Let $D = k_5+k_6+2k_5a_1(1-\beta)S_m(1-E_2).$ We have

Thus,

where $D_1 = (a_1(1-\beta)+S_m(1-E_1))(a_1(1-\beta)+S_m(1-E_2)).$

Now, we consider a special scenario, namely, $E_1 = 0 \; \text{and}\; E_2 = 2E_m.$ Then we can obtain

where $D_1^{s} = (a_1(1-\beta)+S_m)(a_1(1-\beta)+S_m(1-2E_m)).$

Next, we mainly consider the favorable conditions for periodic acaricides spraying in this special case over constant spraying by comparing $R_{ps}^{0s}$ and $R_{cs}^{0}.$ Note that the average acaricidal effect during one year are both $E_m$ for the constant spraying and the special periodic spraying.

Let $\Psi(E_m)$ be the ratio between $R_{ps}^{0s}$ and $R_{cs}^0,$ namely,

which represents the difference degree between $R_{ps}^{0s}$ and $R_{cs}^0.$ It is clear that when $\Psi(E_m) < 1(>1),$ periodic spraying (constant spraying) is more effective as a controlling strategy. Thus, based on $(5.3)$ and $(5.11),$ we define

where $E_m \in (0, \frac{1}{2}),$ and $\Theta(E_m) > 0(<0)$ means $\Psi(E_m) > 1(<1).$ Note that $2S_mE_m < a_1(1-\beta)+S_m$ and $\Theta(0) = 0.$ Let $\psi(E_m) = (a_1(1-\beta)+S_m)(a_1(1-\beta)+S_m(1-2E_m)),$ and then

which is equivalent to

which implies that $\Theta(E_m)$ is a monotonic decreasing function with respect to $E_m$ at $(0, \frac{1}{2}).$ That is to say, $\Theta(E_m) < 0,$ which is equivalent to $\Psi(E_m) < 1.$ Therefore, we can conclude that periodic spraying than the constant one is more beneficial for controlling the amount of tick populations.

6.   Numerical simulations

We use some numerical simulations to illustrate the theoretical results. We list all parameters in Table 1.

For model $(2.1)$, we set diapause rate $d = 0.45$. When birth rate of tick population $b = 30,$ the net reproduction number is $R_0 = 0.9504.$ The global asymptotic stability of the tick-free fixed point $E_0(0, 0, 0)$ is verified and the tick population gradually die out. Figure 2 shows the solutions of model $(2.1)$ with four different initial values. When $R_0 = 1,$ we can compute that $b = 31.5641.$ The eigenvalues of the tick-free fixed point $E_0$ has an eigenvalues $\lambda = 1.$ According to Theorem 2, the tick-free fixed point $E_0$ is a transcritical bifurcation point. Figure 3 shows that the transcritical bifurcation diagram of normal immature tick populations near the fixed point $E_0$ when $b$ varies in the range of $[0, 40].$ From Figure 3, we can obtain that the fixed point $E_0$ is stable when $b < 31.5641$ and is unstable when $b > 31.5641,$ and the stable coexistence fixed point $E^{*}$ emerges when $b > 31.5641.$ Increasing the birth rate to $b = 50,$ the net reproduction number becomes $R_0 = 1.5841.$ Model $(2.1)$ has a unique coexistence fixed point $E^*(173.1298, 23.4383, 40.8961),$ which is locally asymptotically stable(shown in Figure 4).

Now, we explore the acaricidal effect of spraying acaricides. First, we consider the acaricidal effect with constant spraying, namely, $E_m = 0.73.$ When $b = 30, R_{cs}^0 = 0.9722 < 1,$ namely, the tick-free fixed point $E_{cs}(0, 0)$ is globally asymptotically stable (shown in Figure 5 (a) and (b)). When $b = 50,$ we can obtain $R_{cs}^0 = 1.6204 > 1,$ a unique coexistence fixed point $E_{cs}^*(204.1122, 47.6981)$ emerges and its global asymptotic stability is presented in Figure 5(c) and (d).

Moreover, we consider the 2-periodic spraying and fix $E_1 = 0.66\; \text{and}\; E_2 = 0.8.$ When $b = 30, R_{ps}^0 = 0.9884 < 1.$ The tick-free fixed point is globally asymptotically stable and both immature and mature tick populations gradually die out(see Figure 6(a) and (b)); when $b = 50, R_{ps}^0 = 1.6474 > 1,$ a unique and globally asymptotically stable positive 2-periodic solution $\{(246.1702, 45.7523), (166.2203, 49.8790)\},$ we plot Figure 6(c) and (d) to show this dynamic process.

In order to compare the advantages of periodic spraying with constant spraying, we plot Figure 7 to present the net reproduction numbers $R_{cs}^0$ and $R_{ps}^{0s}$ varying with acaricidal effet $E_m\in (0, 0.5).$ Figure 7 shows that $R_{ps}^{0s}$ is always less than $R_{cs}^0.$ It is clear that with the same average acaricidal effect, periodic spraying is more effective than constant spraying for controlling the number of tick populations. Moreover, we numerically simulate tick population dynamics with 4-periodic and 6-periodic acaricidal effect as all parameters remain constant in Figure 8. Table 2 lists the periodic acaricidal effect values, which guarantee the average acricidal effect is $0.73.$ We find that multiple periodic acaricidal effect has a positive impact on the persistence of tick population.

Finally, we compare model $(2.1)$ and model $(5.7)$ respectively when the acaricidal effect is a constant and a 2-periodic function. When the acaricidal effect is constant, we can find that mature tick population in model $(5.7)$ is always larger than that in model $(2.1)$ in Figure 9(a). When the acaricidal effect is a 2-periodic function, the average number of mature tick population in model $(5.7)$ is also always greater than that in model $(2.1),$ which is shown in Figure 9(b). Therefore, diapause can slow down the development of tick population under the same acaricidal effect.

7.   Discussion

This paper studies a two-stage tick population model based on a difference equation system. Diapause is introduced into the model as an essential ecological process in the growth of tick populations. We investigate the properties of model $(2.1)$ by computing the net reproduction number and analyzing the asymptotic stability of fixed points. We suppose that diapause only occurs in immature ticks during unfavorable environment conditions to ensure their survival. Once immature ticks terminate diapause, they will become active and develop into the next stage. Actually, mature ticks may undergo diapause as well, immature ticks moulted and evolved into adults in autumn after feeding in spring and summer, and the mortality rate of ticks increased in autumn due to environmental effects. In order to resist harsh environment, adult ticks undergo behavioral diapause in autumn until they quit diapause and begin feeding in winter of the same year. Behavioural diapause is also thought to prevent ticks that have moulted from spring-fed ticks from feeding in the autumn of the same year ^[14]. In this paper, we ignore diapause in mature ticks for the model simplicity. It is well known that the diapause of ticks can be divided into developmental diapause (temporary suspension of engorged tick development) and behavioral diapause (interruption of host-seeking activity of unfed ticks). Based on our model, if the development of ticks is further refined, two types of diapause are considered at each stage, which will further help us understand the biological dynamics of tick populations and pathogen transmission.

Furthermore, the birth, death and diapause of tick populations are closely related to seasonal changes in environmental conditions such as photoperiod, temperature and humidity, and if the environment changes suddenly and is not conducive to the growth of the ticks, the ticks will be induced into diapause again before the ticks exit diapause and enter the next stage. Thus, it would be more reasonable to incorporate these seasonal variations in the model.The modeling idea of seasonal variation needs to be further extended in the future work, which brings challenges to the theoretical analysis of the model.

At present, acaricides spraying are the main method to control ticks. This paper provides some constructive suggestions for selecting the appropriate control strategy of acaricides spraying. However, the resistance of ticks to acaricides caused by frequent spraying of the same acaricides will make some acaricides ineffective, resulting in potential outbreak of tick populations. The large-scale use of synthetic acaricides will lead to environmental pollution, and even the toxicity in them will harm human and animals^[37]. Optimizing control methods to minimize the resistance of ticks to acaricides and the harm of toxicity will be a challenge in modeling in the future. Moreover, impulsive systems is a better choice for considering acaricide to control tick population. In recent years, some scholars^[39,40] have conducted in-depth research on nonlinear and delayed impulsive systems and obtained some innovative results, which provides great help for applying these theories to our future application research, and the model may present more complex dynamic behavior, such as bifurcation and chaotic phenomena.

Acknowledgments

This work was supported by the National Natural Science Foundation of China (No.12171074).

Conflict of interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Appendix

where $o(|(I(t), I_d(t), M(t), B(t))|^2)$ means terms of order greater than 2 in the combination of $(I(t), I_d(t), M(t), B(t)).$

where $o(|(u(t), v(t), w(t), \mu(t))|^2)$ means terms of order greater than 2 in the combination of $(u(t), v(t), w(t), \mu(t))$ and

It is worth notice that $c_{11} \neq 0.$