Implementation of circular economy principles in management of end-of-life tyres in a developing country (Nigeria)

1.   Introduction

Our focus here is the population dynamics of ticks, such as Ixodes ricinus and Ixodes scapularis, which are responsible for transmitting tick-borne diseases including Lyme disease and tick-borne encephalitis. The lifecycle for ticks consists of four main stages: eggs, larvae, nymphs and adults and the stage-to-stage development from larvae to nymphs and from nymphs to adults occurs after ticks blood feed on hosts. There are different substages for each post-egg stage: the questing stage, the feeding stage and the engorged stage. In the questing stage, ticks are looking for a host to climb on and to feed on. In the feeding stage, they are feeding on a host before detaching and molting on the ground. There is a final stage for female ticks which is the egg-laying stage where ticks lay eggs before dying soon after.

Tick population dynamics are highly dependent on the habitat's local microclimatic conditions ^[1]. Host abundance is key to sustaining the tick population since ticks have to take blood meals from hosts in order to develop. Ticks at different stages of their life cycles have different host preferences: immature ticks mostly feed on small mammals and rodents while adult ticks tend to have their blood meals from large mammal hosts such as deer ^[2]. The abundance of suitable large mammals in the habitat is therefore the key for adult ticks to complete their life cycles and for the reproduction rate of female ticks ^[3]. Importantly, these large mammal hosts also provide the mobility of adult ticks as feeding adult ticks may be carried by the hosts and may drop off in different locations.

Climatic conditions affect the tick population dynamics. Temperature, for example, has been seen to play a crucial role in tick persistence ^[4]. Low temperatures may lead ticks to a quiescence period, a state of torpidity in which ticks arrest their development and slow down their progression ^[5]. It is also well known that ticks cannot survive at really cold temperatures, so when temperature rises due to climate change, the geographic range of ticks expands quickly leading to the observed trend of northward spread of ticks in Canada ^[6]. Other factors influencing tick development include humidity and daylight exposure. Ticks tend to survive easily in humid environments such as those with thick mat layers in the soil ^[7]. Several experimental studies show that ticks are very vulnerable to dessiccation, therefore they remain questing for a limited time before replenishing their reserves ^[8]. It is also well known that ticks are sensitive to changes of daylight and tend to undergo long diapausing periods in habitats where the amount of light is below an optimal level ^[9].

It is therefore plausible that there are, in a given region, different geographical locations with different tick concentrations corresponding to different combinations of local environmental conditions ^[10]. A habitat configuration analysis has been conducted to understand and predict tick spread in a small-scale landscape, and to understand what are the potentially favourable areas for tick development and how connectivity between different patches impacts the tick population distribution ^[11]. In this study, hosts play a key role since ticks themselves have very limited mobility, and they can move between different patches only when they are attached to hosts during their blood meals. Geographical areas can be classified, according to certain micro-habitat and climatic conditions, including those aforementioned, into being favourable or unfavourable for tick growth and spread ^[12]. Certain woodland areas are both less exposed to daylight and contain a dense shrub layer that can keep the humidity level high, these areas are highly favourable for tick growth ^[1]. On the other hand, grasslands do not have any protection from sunlight and can be drier comparing with woodlands, and they are shown to be less favourable for tick development ^[13]. Several studies have found a higher density of ticks in woodlands than in grasslands, for example in Spain ^[14] and Sweden ^[15] for Ixodes ricinus, and in the United States ^[16] for Ixodes scapularis. All the ecological factors mentioned above have been shown to influence tick development by modifying not only their survival rates but also their diapause probabilities ^[17]. This means that ticks in general take less time to develop in favourable areas with respect to unfavourable areas, and have a higher survival rate.

Therefore, human interventions can alter the tick population dynamics. The two main types of control are habitat-related and host-related. The former are used to decrease the favourability of a specific environment for ticks and include the removal of the leaf litter layer, which is important for ticks to avoid dessiccation ^[18], and controlled burns in some tick-rich areas ^[19]. The latter are useful to modify tick migration between patches by fencing deer, for example ^[20]. Both control strategies have a significant impact on the parameters of any tick population dynamics model including the one we are developing here. One of our objectives in this study is to evaluate the impact of these changes on tick population dynamics. Our previous study ^[21] shows that not every control intervention can achieve its intended goals.

Spatial models involving ticks have been extensively used in literature to show how tick population grow and tick-borne diseases spread. Some of these formulations use partial differential equations and aim to study what is the velocity of epidemics spread using traveling waves ^[22,23]. In our approach, we consider a patchy environment instead of a continuous spatial model and focus on the importance of both patch-dependent and host mobility parameters. We will use delay differential equations to capture the physiological structure of ticks. We note that a few delay differential equation compartmental models have been developed. Some incorporated just a single development delay ^[24,25,26]. Some others considered the possibility of ticks to undergo diapause which would lead to an additional delay during development ^[27,28]. The multi-patch approach has also been developed by references ^[29,30] in a multi-species epidemic model and recently by reference ^[31] using an ordinary differential equation system where the focus was mainly on the effect of cofeeding and host movement for disease spread.

The paper is organized as follows. We first introduce the model and show some of its key properties. We then calculate how patch-dependent survival probabilities and migration parameters affect the isolated and interconnected tick reproduction numbers and show their relevance in the study of the equilibria and stability of the model. Finally, we discuss the effect of tick reproduction numbers by showing some simulations and describe the implications of these results.

2.   The patchy model spatially stratified by developmental delay and connected by host mobility

We consider a simplified habitat configuration with patch stratified by tick development delays (regular vs diapause) and connected by host mobility. We assume the region consists of two patches, which are distinguished by the length of life cycle of inhibiting ticks, and connected by large size mammals which provide blood meals to feeding adult ticks. We consider the case where the density of relevant hosts in each patch remains constant. Large size mammals move between the two patches and therefore engorged adult ticks can drop off to both patches due to the host mobility. We ignore the mobility of hosts for larval and nymphal ticks, so in our patchy model, only feeding adult ticks and egg-laying adult ticks are explicitly incorporated while ticks in other stages can be calculated from the production rate and the survival probability. To quantify the movements of hosts for feeding adult ticks between two patches, we assume a portion $\alpha _{ij}$ of feeding adult ticks in patch $i$ can drop off to patch $j$ to become egg-laying adults in patch $j$ and we consider this movement to be instantaneous. We consider the case where ticks grow up from the eggs in a given patch remain in the patch until they reach the stage of feeding adults, in other words, hosts for larval and nymphal ticks can move within the patch but not to the other patch. We allow the two patches to be distinct in terms of the developmental delays from eggs to feeding adults, one with normal developmental delay $\tau _1$ and another with diapause developmental delay $\tau _2$. Since diapause would result in a longer overall development time for ticks, we assume $\tau_2 > \tau_1$. The probability of survival from eggs to feeding adults $\rho_i$ depends exclusively on the patchy environment, although this can be relaxed to allowing delay-dependent probability. The birth rate is, as normal, assumed to be dependent on the egg-laying adult density in the same patch by the Ricker reproduction function $f(x) = pxe^{-qx}$ ^[32] with $p$ and $q$ being positive constants. The transition time from feeding adult ticks to engorged ticks is relatively short, so the probability for the engorged adult ticks to drop off to the same patch where the feeding ticks come from is relatively large, and hence $\alpha_{12} + \alpha_{21} < 1$. In our model, feeding adults either die out with the death rate $\mu$ or develop into egg-laying adults with the development rate $\gamma$; and feeding adults advance to the egg-laying stage with rate $\gamma_F = \gamma \theta$, where $\theta$ represents the survival probability of female adult ticks in the engorged stage. Egg-laying adults then die out with the death rate $\delta$. The parameters are incorporated in Figure 1 and have been summarized in Table 1.

With the above assumptions and notations, we can now formulate the patchy model with multiple delays as follows:

In the above formulation, we use $F_i(t)$ to denote the number of feeding adults and $L_i(t)$ for the number of egg-laying adults in patch $i$.

3.   Preliminary qualitative properties of the model

We first show that solutions to (2.1) are non-negative (given non-negative initial conditions) and bounded. Define the phase space

with norm

Proposition 1. Solutions to (2.1) (given non-negative initial conditions) are non-negative and uniformly bounded.

Proof. For a given initial data $\phi \in X^+$, we obtain a unique solution to (2.1) $x^{\phi}(t)$ for positive times (^[42], Theorem 2.3), generating a semiflow on $X^+$. We can easily show that the solution to (2.1) with initial data $\phi \in C([-\tau_2, 0], \mathbb{R}^4_+)$ remains non-negative. In what follows, we show that solutions of the model with non-negative initial data are attracted to a bounded and positively invariant subset in $X^+$. Let

where

We note that the Ricker function $f(x) = pxe^{-qx}$ is bounded for $x \geq 0$ and has its maximum at $x = \frac{1}{q}$. Let

Therefore,

from which and with $F_1(0) = F_1^0$ and $F_2(0) = F_2^0$ it follows that

In particular, if $F_i^0 \leq F_i^\infty$, then $F_i(t) \leq F_i^\infty, \; \forall t \geq 0, \ {\rm{ for }}\ i = 1, 2$. Also note that

from which and with $L_1(0) = L_1^0$ and $L_2(0) = L_2^0$ it follows that

Consequently, if $L_i^0 \leq L_i^\infty$, then $L_i(t) \leq L_i^\infty \; \forall t \geq 0, \ {\rm{ for }}\ i = 1, 2$. The above argument implies that all solutions of the model system remain bounded for all $t\ge 0$ and solutions are in fact ultimately uniformly bounded since

Let

Then we have shown that $X^+_\Gamma$ is a positively invariant set in $X^+$ which attracts solutions of (2.1) with initial data in $X^+$.

4.   Tick reproduction number

We start with a special case where two patches are isolated from each other, and we have for $i = 1, 2$ the coupled system

We compute the tick basic reproduction number $R_{0, i}$ for a single patch $i$ which allows us to study the average number of female ticks that are born by a single female tick in this patch. The procedure is similar to calculating $R_0$ for epidemics where we consider linearization at the tick-free equilibrium. The existence of a positive feedback $f'(0) = p > 0$ guarantees that monotone dynamical theory can be applied. Therefore, the stability of the trivial solution of (4.1) is equivalent to that of the linear ordinary differential equation system associated (where $\tau_i$ = 0) ^[33].

We first linearize system (4.1) at the trivial equilibrium, where $\tau_1 = \tau_2 = 0$, to get

with the transmission matrix $T$ and the transition matrix $\Sigma$

where $p = f'(0)$. To use the next generation matrix approach, we note that the inverse of $\Sigma$ and the next generation matrix $K = -T \Sigma^{-1}$ are given by

Therefore the spectral radius

gives the so-called basic reproduction number. As shall be shown, this basic reproduction number decides whether ticks can persist in patch $i$, namely, $R_{0, i} > 1$ implies persistence and $R_{0, i} < 1$ implies extinction of ticks in the patch. In what follows, we focus on the case in which $R_{0, 1} > 1 > R_{0, 2}$. So when the two patches are isolated from each other, ticks can persist in patch 1 but will be extinct in patch 2.

We now consider the case when two patches are connected by mobility of feeding adults ($\alpha_{12}, \, \alpha_{21} \neq 0$). Similarly to the isolated case, the stability of the trivial solution of (2.1) is equivalent to that of the linear ordinary differential equation system associated (where $\tau_1 = \tau_2$ = 0) ^[33]. We linearize model (2.1) at the trivial equilibrium to get a linear system of ordinary differential equations for $X = (F_1, F_2, L_1, L_2)^T$, with $\tau_1 = \tau_2 = 0$

We now use the next generation approach to introduce the so-called basic reproduction number. The transmission matrix $T$ and the transition matrix $\Sigma$ are given by

respectively.

The inverse of $\Sigma$ is given by

so the next generation matrix $K = -T \Sigma^{-1}$ is given by

We define the tick reproduction number $R_{0, c} = \rho(K_{12})$, where $\rho$ denotes the spectral radius and

Therefore the characteristic equation of $K_{12}$ is:

and can be rewritten as

where

Note that both eigenvalues are real since

Therefore the tick reproduction number is

There are two interesting special cases of semi-connectedness of two patches.

4.1.   Escalating up

This is the case when hosts for feeding adults move only from patch 2 (tick-unfavorable patch) to patch 1 (tick-favorable patch), so $\alpha_{12} = 0$. In this case, we have

and

4.2.   Cascading down

This is the case when hosts for feeding adults move only from patch 1 (tick-favorable patch) to patch 2 (tick-unfavorable patch), so $\alpha_{21} = 0$. In this case, we have

and

5.   Model equilibrium analyses

A nontrivial equilibrium $(F_1^*, F_2^*, L_1^*, L_2^*)$ of (2.1) is given by

This can be rewritten as:

5.1.   Isolated patches

Without host mobility $(\alpha_{12} = \alpha_{21} = 0)$, (5.1) is given by

Using the basic reproduction number, we obtain

Noting $\frac{f(x)}{x} = p e^{- qx}$, we get

Recall that $\rho_1 \geq \rho_2, \, R_{0, 1} \geq R_{0, 2}$, we conclude that if $R_{0, 1} < 1$, then the unique equilibrium is $E_0 = (0, 0, 0, 0)$; if $R_{0, 2} < 1 < R_{0, 1}$, then the model has a nontrivial equilibrium $E_1 = (F_1^*, 0, L_1^*, 0)$; if $R_{0, 2} > 1$, then the model has three nontrivial equilibria: $E_1 = (F_1^*, 0, L_1^*, 0)$, $E_2 = (0, F_2^*, 0, L_2^*)$, and the coexistence equilibrium $E_C = (F_1^*, F_2^*, L_1^*, L_2^*)$. As stated before, we focus on the second case. In addition to the threshold $R_{0, i} = 1$, there is another threshold $R_{0, i} = e$ at which the system changes its feedback nature from positive ($R_{0, i} < e$) to negative ($R_{0, i} > e$). At this point, the resulting equilibrium $L_i^* = \frac{1}{q}$ maximizes the birth function $f(L_i)$ and separates the cases for which $f'(L_i^*) > 0$ (for $R_{0, i} < e$) and $f'(L_i^*) < 0$ (for $R_{0, i} > e$).

5.2.   Semi-connected patches

In the escalating up case when ticks move only from an unfavourable to a favourable environment (i.e., $\alpha_{12} = 0$), there exists a unique non-trivial equilibrium $(F_1^*, 0, L_1^*, 0)$ where $(F_1^*, L_1^*)$ is identical to the isolated case in the scenario considered $(R_{0, 2} < 1 < R_{0, 1})$. In the cascading down case (i.e., $\alpha_{21} = 0$), we have either only the trivial equilibrium or a coexistence equilibrium. Namely, if $\alpha_{12} > \frac{R_{0, 1} - 1}{R_{0, 1}}$, we have the equilibrium $E_0 = (0, 0, 0, 0)$, but when $\alpha_{12} < \frac{R_{0, 1} - 1}{R_{0, 1}}$, we have a coexistence equilibrium $E_C = (F_1^*, F_2^*, L_1^*, L_2^*)$ that will be specified below.

Indeed, from the non-trivial equilibrium equations in (5.1), we get

From the first equation of (5.2) and using $\frac{f(x)}{x} = p e^{- qx}$, we find that $L_1^* = q^{-1} \ln((1-\alpha_{12})R_{0, 1}) > 0$ only if $\alpha_{12} < \frac{R_{0, 1} - 1}{R_{0, 1}}$. In this case we rewrite the second equation of (5.2) as

Noting that $R_{0, 2} = \frac{\rho_2 p \gamma_F}{\delta (\gamma + \mu)}$ and defining $\zeta^* : = \frac{\alpha_{12}}{1-\alpha_{12}} L_1^* \frac{p}{R_{0, 2}}$, we have

and $\zeta^* > 0$ guarantees that there is always a nontrivial solution of this equation.

5.3.   Interconnected patches

We now consider the fully connected case. First of all, using the equations in (5.1), we find $F_1^*$ and $F_2^*$ as linear combinations of $L_1^*$ and $L_2^*$ whenever $\alpha_{12} \neq 0$, $\alpha_{21} \neq 0$ and $\alpha_{12} + \alpha_{21} < 1$. These are given by

where $a_{11} = \frac{\delta (1-\alpha_{21})}{\gamma_F (1 - \alpha_{12} - \alpha_{21})}$, $a_{12} = - \frac{\delta \alpha_{21}}{\gamma_F (1 - \alpha_{12} - \alpha_{21})}$, $a_{21} = - \frac{\delta \alpha_{12}}{\gamma_F (1 - \alpha_{12} - \alpha_{21})}$, $a_{22} = \frac{\delta (1 - \alpha_{12})}{\gamma_F (1 - \alpha_{12} - \alpha_{21})}$.

Note that $\alpha_{12} + \alpha_{21} < 1 \ {\rm{ if\ and\ only\ if }}\ a_{11}a_{22} > a_{12}a_{21}$.

5.3.1.   Geometric approach

We now develop a geometric approach to look at (5.3) as we change $\alpha_{12}$ and $\alpha_{21}$ subject to the constraint $\alpha_{12} + \alpha_{21} < 1$. For notation simplicity, let

Then (5.3) becomes

So we substitute the feeding adult equilibria in the two patches $F_1^*$ and $F_2^*$ in (5.3) with the first two equations of (5.1) and get

We want to explore graphically the behaviour of $L_1^*$ with respect to $L_2^*$, to study conditions for coexistence equilibrium, which would be the intersection of the functions defined below. Taking $L_2^*$ as a function of $L_1^*$ yields

and $L_1^*$ as a function of $L_2^*$ yields

In particular we have $F(x) - \frac{1 - \alpha_{21}}{\alpha_{21}} x \rightarrow 0$ as $x \rightarrow \infty$ and $G(x) - \frac{1 - \alpha_{12}}{\alpha_{12}} x \rightarrow 0$ as $x \rightarrow \infty$. In order for $F$ and $G$ to be plotted on a $x = L_1^*, y = L_2^*$ plot, we need to reflect $G$ about the line $y = x$. Note also that

Therefore $F'$ and $G'$ are increasing functions since $f'$ is decreasing, with

Recall that we are interested in the case when $R_{0, 1} > 1 > R_{0, 2}$. This means that $F'(0)$ is relatively small while $G'(0)$ is relatively large. Since $0 < R_{0, 2} < 1$, we see that $G(x)$ is always non-negative for $x > 0$, since $G(0) = 0$, $G'(0) > 0$ and $G'(x)$ is an increasing function for $x \geq 0$. Therefore $G^{-1}$ is well defined for $x \geq 0$, and $(G^{-1})'(x)$ is a decreasing function for $x \geq 0$, and $G^{-1}(x) - \frac{\alpha_{12}}{1 - \alpha_{12}} x \rightarrow 0$ as $x \rightarrow \infty$. We are interested in studying the possible intersections between $F(x)$ and $G^{-1}(x)$. We see that $G^{-1}(x) < F(x)$ for large $x$ since the following inequality holds under the condition $\alpha_{12} + \alpha_{21} < 1$:

The only case there can be one and only one non-trivial intersection between $F(x)$ and $G^{-1}(x)$ occurs if

Therefore, we introduce the threshold value $T_{coex} : = F'(0) G'(0)$

In particular, the existence of a non-trivial equilibrium occurs if $T_{coex} < 1$.

Recall that we have computed another threshold value $R_{0, c}$ in (4.3), which is also a threshold to determine if a coexistence equilibrium exists. The next result shows the equivalence of $R_{0, c}$ and $T_{coex}$ in terms of the coexistence equilibrium.

Theorem 5.1.

Proof. From the expression of $R_{0, c}$, we note that

We consider two cases: the case $b > 2$, and hence $\sqrt{b^2 - 4c} \geq 0 > 2 - b$ and the case $b \leq 2$. In this case $R_{0, c} > 1\ {\rm{ if\ and\ only\ if }}\ b > c+1$.

In a similar, way we study when $T_{coex} < 1$, which is equivalent to

which can be rewritten as

By dividing both sides by $(1-\alpha_{12}-\alpha_{21}) > 0$, we have that $b > c+1$. Note that this condition is equivalent to that for $R_{0, c} > 1$ in case $b \leq 2$.

It remains to show that $b > 2$ implies $T_{coex} < 1$. We know $R_{0, c} \in \mathbb{R}^+$, therefore $b^2 \geq 4c$. But $b > 2$ means $c \leq 1$; so $b > c+1$ holds for $b > 2$ and $c \leq 1$. We can also use a similar argument to show $R_{0, c} = 1\ {\rm{ if\ and\ only\ if }}\ T_{coex} = 1$.

So we have the following result.

Corollary 1. If $R_{0, c} < 1 \, (\mathit{{i.e., }}\ T_{coex} > 1)$, then the model has only the trivial equilibrium $(0, 0, 0, 0)$. If $R_{0, c} > 1 \, (\mathit{{i.e., }}\ T_{coex} < 1)$, then the model has a nontrivial equilibrium $(F_1^*, F_2^*, L_1^*, L_2^*)$ with each component positive.

5.3.2.   Asymptotic expansion analyses

We have seen that in the interconnected case a closed form solution for the coexistence equilibrium is difficult to obtain. In this subsection, we study the coexistence equilibrium using a perturbation analysis when host mobility is small. The equilibria in the case when $\alpha_{12} = \epsilon \alpha_{12}^0$ and $\alpha_{21} = \epsilon \alpha_{21}^0$ with $\epsilon < < 1$ can be computed by solving

Consider the asymptotic expansion of the equilibria as:

where $(f_1, f_2, l_1, l_2) = (\frac{\delta}{\gamma_F q} \ln(R_{0, 1}), 0, \frac{1}{q} \ln(R_{0, 1}), 0)$ is the equilibrium in the isolated patches case where $R_{0, 1} > 1 > R_{0, 2}$. Making the appropriate substitutions, (5.5) becomes

By substituting the values in (5.4), we have the following equations:

Ignoring $o(\epsilon)$, this is a linear system of four equations and four unknowns $l_{1, 2}, \, l_{2, 2}, \, f_{1, 2}, \, f_{2, 2}$ which yields the solution

We therefore obtain the final asymptotic expansion of the coexistence equilibrium in (5.6) as

We observe that, first of all, the solution $(L_1^*, L_2^*, F_1^*, F_2^*)$ is always positive since $R_{0, 1} > 1 > R_{0, 2}$ and $\epsilon < < 1$. Secondly, a small host mobility causes the decrease of egg-laying adults in the favorable environment and an increase in the unfavorable patch. Also, the feeding adult equilibrium in the favorable environment can increase or decrease, depending on the size of $R_{0, 1}$. If $R_{0, 1} > e$, then $F_1^*$ increases, otherwise if $R_{0, 1} < e$, then $F_1^*$ decreases. Finally, the fact that $F_1^*$ increases with mobility of the host if $R_{0, 1} > e$ is interesting: although mobility of the host decreases $L_1^*$, a large basic reproduction number in the favorable patch amplifies the number of ticks in the stage from eggs to feeding adults to compensate for the loss of $L_1^*$.

6.   Stability

6.1.   Linearization

Linearization at the trivial equilibrium $(0, 0, 0, 0)$ yields:

The characteristic equation for $\lambda \in \mathbb{C}$ is given by $\det(B) = 0$ where $B = B_0 + B_1 e^{-\lambda \tau_1} + B_2 e^{-\lambda \tau_2} - \lambda I$,

and

Defining $\psi: = (\gamma + \mu + \lambda) (\delta + \lambda)$, we have

where

Similarly, the linearization at the non-trivial equilibrium $(F_1^*, F_2^*, L_1^*, L_2^*)$ is given by

The characteristic equation derives from $\det(C) = 0$, where

In particular, $\det(C) = \psi^2 - c_0 \psi + c_1$ and

6.2.   Isolated patches

With $\alpha_{12} = \alpha_{21} = 0$, we have

Therefore $\det(B) = 0\ {\rm{ if\ and\ only\ if }}\ (\psi - \gamma_F \rho_1 p e^{-\lambda \tau_1}) (\psi - \gamma_F \rho_2 p e^{-\lambda \tau_2}) = 0$. Here, each factor corresponds to the characteristic equation of an isolated patch linearized at the trivial equilibrium. Similarly we have $\det(C) = 0$ written as

Proposition 2. The trivial equilibrium $(0, 0)$ of (4.1) is a global attractor if $R_{0, i} < 1$ and unstable if $R_{0, i} > 1$.

Proof. The linearized system of patch $i$ at the trivial equilibrium is given by

This is a delay differential system with an irreducible and cooperative delayed feedback. By Corollary 5.2 of Smith ^[33], the stability of the above system is the same as that of the corresponding ordinary differential equation model (letting $\tau_i = 0$) from which the conclusion follows.

We can also consider the stability of the non-trivial equilibrium $(F_i^*, L_i^*) = \left(\frac{\delta \ln(R_{0, i})}{\gamma_F q}, \frac{\ln(R_{0, i})}{q}\right)$, when exists, by considering the linearization

where

It has been shown that $1 < R_{0, i} < e \Rightarrow f'(L_i^*) > 0$; and $R_{0, i} > e \Rightarrow f'(L_i^*) < 0$. Let

The characteristic equation is $\det B = 0$, namely,

and can be rewritten as

Note that if $1 < R_{0, 1} < e$, the system is cooperative and the stability of non-trivial equilibrium of (6.3) is same for all $\tau_i$ so we can impose $\tau_i = 0$ and study the ODE model instead ^[33]. We show that local stability of the non-trivial equilibrium holds also when $e < R_{0, 1} < e^2$.

Proposition 3. The non-trivial equilibrium $(F_i^*, L_i^*)$ of (4.1) is locally asymptotically stable for all $\tau_i \geq 0$ if $1 < R_{0, i} < e^2$.

Proof. We want to show that there are no solutions to the characteristic equation (6.6) with positive real part. Suppose by contradiction there exists a root of (6.6) $\lambda = x + iy$ with $x \geq 0$. So the following equality holds:

We also know

While the right hand side with $1 < R_{0, i} < e^2$ satisfies

a contradiction.

Proposition 4. Every non-trivial solution of (4.1) converges to $(F_i^*, L_i^*)$ as $t \rightarrow \infty$ if $1 < R_{0, i} < e$.

Proof. We have shown that $\Gamma$ defined in (3.1) is a positively invariant set of (2.1) and the $\omega$-limit of the solutions is in $\Gamma$. Let $\Gamma_{iso} : = \{(F_i, L_i) \in \mathbb{R}^2_+ : F_i \leq F_i^{\infty}, L_i \leq L_i^{\infty}\}$, where $F_i^{\infty} = \frac{\rho_i p}{q e (\gamma + \mu)}, L_i^{\infty} = \frac{\gamma_F p \rho_i}{q e \delta (\gamma + \mu)}$. Consider the Jacobian of (6.3) for $\tau_i = 0$:

We see that $\Gamma_{iso}$ is a positively invariant set of (4.1) containing $(F_i^*, L_i^*)$ and the $\omega$-limit of its solutions is contained in $\Gamma_{iso}$, using the argument similar to that for Proposition 1. We also note that (6.3) is a cooperative system (for $\tau_i = 0$), since $j_{1, 2}, \, j_{2, 1} \geq 0$ as $f'(L_i^*) > 0$ if $1 < R_{0, i} < e$. Finally, $J$ is irreducible since $j_{1, 2}, \, j_{2, 1} \neq 0$. Therefore, using the monotone dynamical systems theory ^[33], we conclude that $(F_i^*, L_i^*)$ is globally attractive.

We remark that it is possible to extend global attractivity of the non-trivial equilibrium also for $e \leq R_{0, i} < e^2$ using exponential ordering ^[33].

Proposition 5. Assume $R_{0, i} > e^2$. The equilibrium $(F_i^*, L_i^*)$ of (4.1) is locally asymptotically stable if $\tau_i < \tau^*$, and is unstable for $\tau_i > \tau^*$, where $\tau^* = \omega^{-1} \arctan\left(\frac{\omega (\delta + \gamma + \mu)}{\omega^2 - \delta (\gamma + \mu)}\right)$.

Proof. We start with the asymptotic stability of the non-trivial equilibrium when $\tau_i = 0$. The characteristic equation becomes $(\gamma + \mu + \lambda)(\delta + \lambda) = \rho_1 \gamma_F f'(L_i^*)$, which can be rewritten as

It is easy to check that both zeros have negative real parts, so $(F_i^*, L_i^*)$ is asymptotically stable in the absence of delay.

We expect the presence a Hopf bifurcation as the delay grows. In order to apply the Hopf bifurcation theorem for DDE systems ^[34], suppose there exists a purely imaginary root $\lambda^* = \pm i \omega$. Then $(\gamma + \mu + i \omega)(\delta + i \omega) = \rho_1 \gamma_F f'(L_i^*) e^{- i \omega \tau_1^*}$. We separate the real and imaginary parts of this system and yield

The Hopf bifurcation point $\tau_i^*$ is the smallest $\tau_i$ that satisfies (6.7). Summing up the square of first line to the square of the second line yields

Letting $\zeta = \omega^2$, we obtain

That is, using (6.4), we have

Note that if $R_{0, i} > e^2$, there will always be a unique positive solution since

The two real solutions of $\zeta$ are $\zeta_{1, 2} = \frac{1}{2} (2 \delta (\gamma + \mu) - (\delta + \gamma + \mu) \pm \sqrt{\Delta})$, and the only positive solution is

Therefore $\omega = \pm \sqrt{\zeta_1}$. To find $\tau_i^*$, use (6.7) to get

therefore

Consider

Since $\lambda^*$ is a root of (6.5) for $\tau_i = \tau_i^*$, $h(\lambda^*, 0) = 0$. We can verify that $h_{\lambda}(\lambda^*, 0) \neq 0$, so the Implicit Function Theorem ensures that there exists $\lambda(\xi) \in C^1$ such that $\lambda(0) = \lambda^* = i\omega$. So proving transversality condition reduces to showing that ${\rm{Re}}(\lambda'(0)) \neq 0$. Using implicit differentiation and noting that $h(\lambda(\xi), \xi) = 0$, we have

The computation of the partial derivatives of $h$ at $(\lambda^*, 0)$ yields

With $\lambda^* = i\omega$ and using the fact that ${\rm{Re}}(z) = \frac{ac+bd}{c^2-d^2}$ for a complex number $z = \frac{a+ib}{c+id}$ where $(c, d)\neq (0, 0)$, and using the Euler's formula, and after a series of computations, we have

So, we need to show that $\tan (\omega \tau_i^*)(\gamma + \mu +\delta) + 2\omega \neq 0$, which by the computed tangent from (6.10), is equivalent to $(\delta + \gamma + \mu)^2 + 2[\omega^2 - \delta(\gamma+\mu)] \neq 0$. This, by isolating $\omega^2$, is equivalent to $\omega^2 \neq \frac{1}{2}(2\delta(\gamma + \mu) - (\delta + \gamma + \mu)^2)$. But $\omega^2 = \zeta_1$ in (6.9) so we want to compare the two quantities and check if $\sqrt{\Delta} \neq (\delta + \gamma + \mu) - (\delta + \gamma + \mu)^2$. Using (6.8) for $R_{0, i} > e^2$ and the fact that $(\delta + \gamma + \mu)^2 \geq 2 \delta(\gamma + \mu)$, we get

completing the proof of the transversality condition.

We remark that the critical value $\tau_i^*$ decreases when $p \rho_i \gamma_F$ increases, therefore survival probabilities and maximal amount of eggs produced influence the Hopf bifurcation. Also $\omega^2$ depends directly on $b: = \delta + \gamma + \mu -2 \delta (\gamma + \mu)$. In this case as $b$ decreases, also $\tau_i^*$ decreases. $(F_1^*, 0, L_1^*, 0)$ is asymptotically stable if $1 < R_{0, 1} < e^2$ for all $\tau_1 > 0$; or $R_{0, 1} > e^2$ and $\tau_1 < \tau^*$; and unstable if $\tau_1 > \tau^*$.

6.3.   Semi-connected patches

We now consider the semi-connected case.

Escalating up: In this case, $\alpha_{12} = 0$, we have both trivial equilibria and a non-trivial equilibrium $(F_1^*, 0, L_1^*, 0)$. The characteristic equation at the non-trivial equilibrium $\det(B) = \psi^2 - c_0 \psi + c_1 = 0$, where

We can factorize $\det(B) = 0$ as $(\psi-a)(\psi-b) = 0$, where $a = \gamma_F \rho_1 f'(L_1^*) e^{- \lambda \tau_1}$ and $b = (1-\alpha_{21}) \gamma_F \rho_2 f'(L_2^*) e^{- \lambda \tau_2}$.

Patch 2 only admits the equilibrium $(0, 0)$ which is asymptotically stable. The stability of the equilibria in Patch 1 instead will be exactly the same as that of the isolated case, therefore we conclude that: trivial equilibrium $(0, 0, 0, 0)$ is unstable (Proposition 2); the non-trivial equilibrium $(F_1^*, 0, L_1^*, 0)$ is asymptotically stable if $1 < R_{0, 1} < e^2$ for all $\tau_1$ (Proposition 4); is asymptotically stable if $R_{0, 1} > e^2$ for $\tau_1 < \tau^*$ and unstable otherwise since it undergoes a Hopf bifurcation at a certain $\tau^*$ derived by (6.7) (Proposition 5).

Cascading down: In this case, $\alpha_{21} = 0$, there is either only the trivial equilibrium if $\alpha_{12} > \frac{R_{0, 1} - 1}{R_{0, 1}}$, or there is also a coexistence equilibrium $(F_1^*, F_2^*, L_1^*, L_2^*)$.

The characteristic equation at the trivial equilibrium is $\det(B) = \psi^2 - c_0 \psi + c_1$ with

Thus, we have $\det(B) = (\psi-a)(\psi-b)$ with $a = (1-\alpha_{12}) \gamma_F \rho_1 f'(L_1^*) e^{- \lambda \tau_1}$ and $b = \gamma_F \rho_2 f'(L_2^*) e^{- \lambda \tau_2}$. Note that $\psi-b = 0$ will results only in eigenvalues with negative real part since $f'(L_2^*) \leq p$ and $R_{0, 2} < 1$, we need to only focus on $\psi-a = 0$. Therefore, we have: the trivial equilibrium $(0, 0, 0, 0)$ is asymptotically stable if $(1 - \alpha_{12}) R_{0, 1} < 1$, and unstable if $(1 - \alpha_{12}) R_{0, 1} < 1$; the coexistence equilibrium $(F_1^*, F_2^*, L_1^*, L_2^*)$ exists if $(1 - \alpha_{12}) R_{0, 1} < 1$ and is asymptotically stable if $1 < (1 - \alpha_{12}) R_{0, 1} < e^2$ for all $\tau_1$; asymptotically stable if $(1 - \alpha_{12}) R_{0, 1} > e^2$ for $\tau_1 < \tau^*$; and unstable otherwise since it undergoes a Hopf bifurcation at $\tau_1 = \tilde{\tau}$ derived below.

Hopf bifurcation occurs when there exists a purely imaginary root $\lambda = i \omega$ to $\psi-a = 0$, which is found by solving the following system:

Using calculations similar to the isolated patch case in Proposition 5, we have

6.4.   Interconnected patches

Stability in the interconnected case presents several complications due to the fact that the equilibrium is not in a closed form and that the characteristic equation cannot be factored out. We analyse the solutions of the characteristic equation

where

and

in terms of the basic reproduction numbers $R_{0, 1}$ and $R_{0, 2}$.

First, we study the solution of this second order equation with respect to $x$, therefore we compute

So

where

and we can rewrite the characteristic equation as

Therefore by replacing $\psi$ with its definition, we have

and the two different terms can be rewritten separately for $i = 1, 2$ as

Therefore the solutions of the characteristic equation (i.e. the eigenvalues of linearized matrix at trivial equilibrium) have to satisfy one of the following 4 equations for $i = 1, 2$:

where

Note that $\psi_i$ is a function of $\lambda$. At this point we prove an analogous result to the isolated case.

Proposition 6. The trivial equilibrium $(0, 0, 0, 0)$ of (2.1) is a global attractor if $R_{0, c} < 1$ and unstable if $R_{0, c} > 1$.

Proof. The linearized system of patch i at the trivial equilibrium is given by (6.1). This is a delay differential system with an irreducible and cooperative delayed feedback. By Corollary 5.2 of Smith ^[33], the stability of the above system is the same as that of the corresponding Ordinary Differential Equation model

From that, the conclusion follows.

Analytical considerations on the coexistence equilibrium cannot be made since there is no closed-form solution. We can study though the stability of the coexistence equilibrium in the asymptotic case where $\epsilon < < 1$

and show that the stability properties of (4.1) are preserved.

Proposition 7. The non-trivial equilibrium $(F_1^*, F_2^*, L_1^*, L_2^*)$ in (5.6) is locally asymptotically stable in (6.13) for all $\tau_1, \, \tau_2 \geq 0$ if $1 < R_{0, c} + o(\epsilon) < e^2$.

Proof. We calculate $R_{0, c}$ for (6.13)

where

Through a series of algebraic computations we have that

We consider the asymptotic expansion for $x \rightarrow 0$ and $a > 0$,

and observe that

So the asymptotic form of $R_{0, c}$ is

We study the characteristic equation

where

and the parameters are

and

As in the isolated case, we can rewrite the equation as

Therefore the solutions with respect to $x$ are

We want to show there are no solutions to the characteristic equation with positive real part:

Suppose by contradiction there exists a root of (6.15) $\lambda = x + iy$ with $x \geq 0$. So the following equality holds for $j = 1, 2$:

We also know

Taking the asymptotic solution from (5.6) and using the fact that $f(x) = pxe^{-qx}$, we see that

For $j = 1$, the right hand side with $1+\epsilon \alpha_{12} < R_{0, i} < e^2(1+\epsilon \alpha_{12})$ satisfies

Note that the condition $1+\epsilon \alpha_{12}^0 < R_{0, i} < e^2(1+\epsilon \alpha_{12}^0)$ is satisfied if $1 < R_{0, c} < e^2$ using (6.14).

A similar conclusion can be inferred for $j = 2$ when $R_{0, 2} < 1$

We thus reach a contradiction, therefore all the roots of the characteristic equation have negative real part and $(F_i^*, L_i^*)$ is locally asymptotically stable for $1 < R_{0, i} < e^2$.

Proposition 8. Every non-trivial solution of (6.13) converges to the equilibrium $(F_1^*, F_2^*, L_1^*, L_2^*)$ in (5.6) for all $(\tau_1, \tau_2) \geq 0$ if $1 < R_{0, c} + o(\epsilon) < e$.

Proof. We have shown that $\Gamma$ defined by (3.1) is a positively invariant set of (2.1) and the $\omega$-limit of the solutions is in $\Gamma$. In a similar way, by multiplying both migration terms by $\epsilon$, we find a set $\Gamma_{asy}$ preserving the same properties for (6.13). Let

where

Consider the Jacobian of (6.13) for $\tau_1 = \tau_2 = 0$:

We see from (6.16) that for small $\epsilon$, $f'(L_2^*) > 0$. We want to understand how $f'(L_1^*)$ is related to $R_{0, c}$ using (5.6). We study the case in which $f'(L_1^*) > 0$, which corresponds to the case in which $L_1^* < \frac{1}{q}$.

From (6.14) we see that $f'(L_1^*) > 0$ if and only if $R_{0, c} + o(\epsilon) < e$. We see that

1. $\Gamma_{asy}$ is a positively invariant set of (6.13) containing $(F_1^*, F_2^*, L_1^*, L_2^*)$ in (5.6) and the $\omega$-limit of its solutions is contained in $\Gamma_{asy}$ - check Proposition 1 for derivation.

2. (6.13) is a cooperative system (for $\tau_1 = \tau_2 = 0$) since $j_{kl} \geq 0$ for $k \neq l$. This is true since $f'(L_1^*), f'(L_2^*) > 0$ if $1 < R_{0, c}+o(\epsilon) < e$.

3. $J$ is irreducible since for every nonempty proper set $I$ of $N = \{1, 2, 3, 4\}$, there is a $k \in I, \; j \in N\setminus I$ such that $j_{kl} \neq 0$ and the digraph is strongly connected.

Therefore, using monotone dynamical theory ^[33], we deduce that (6.13) does not contain any periodic solution for $1 < R_{0, c} < e$, therefore $(F_1^*, F_2^*, L_1^*, L_2^*)$ is attractive.

It is possible to extend global attractivity of the coexistence equilibrium also for $e \leq R_{0, c} < e^2$ using exponential ordering ^[33] and to study Hopf bifurcation of (6.13).

For the general interconnected case (2.1), we observe through simulations that $R_{0, c}$ has a similar behaviour to $R_{0, i}$ and the properties that have been proved by studying the isolated model and asymptotic expansion can be extended to the general model. In particular we observe that

● $1 < R_{0, c} < e^2$ - the coexistence equilibrium is a global attractor of (2.1).

● $R_{0, c} > e^2$ - the coexistence equilibrium is conditionally asymptotically stable on the choices of $\tau_1$ and $\tau_2$.

7.   Simulations

The parameters considered in the following simulation are mainly deriving from literature and are mostly fixed. The normal development delay $\tau_1$ is set to 2 years (730 days) while the diapause development delay $\tau_2$ is set to 3 years (1095 days). The development time from feeding adult to egg-laying is between 10 and 28 days ^[35] so we choose it to be on average two weeks (i.e., $\gamma = 1/14)$. The death rate of feeding adults is set to $\mu = 0.005$ ^[36,37], the survival probability of ticks from feeding to egg-laying stage is $\theta = 0.81$ ^[38] and the exit rate $\delta$ for feeding adults is set to 1 ^[35]. The final fixed parameters are the two constants of the Ricker functions which are $p = 1000$ and $q = 6.2$ ^[38,39]. We used Matlab to provide the simulations using the dde23 algorithm with initial conditions $F_1 = L_1 = 1$ and $F_2 = L_2 = 0.1$ and the biftool package. Note that the choice of initial conditions does not have an impact on the plots in the long run.

Since the aim of this paper is to study how environment and movement affect tick dynamics, survival probabilities and migration coefficients between both patches will be the only parameters that vary throughout the experiments. Note that $\rho_1$ and $\rho_2$ are always chosen such that the condition of favourable and unfavourable environment are satisfied ($R_{0, 1} > 1$ and $R_{0, 2} < 1$).

7.1.   Isolated patches

In the absence of host migration, ticks always die out in patch 2, while in patch 1 they converge to a non-trivial equilibrium in the right plot or to a periodic solution in the left plot of Figure 2. The difference in these two plots is the choice of the survival probability which is larger in the right plot ($\rho_1 = 0.015$) with respect to the former ($\rho_1 = 0.0065$). We have previously shown that the threshold $R_{0, i} = e^2$ is key to study global dynamics of the model. In the right plot, periodic solutions occur since $R_{0, 1} = 11.36 > e^2$ and the 3-year delay is larger than the threshold delay $\tau_1^*$ which can be computed using Hopf bifurcation theory to be $\tau_1^* \sim 28$ days in this specific case.

7.2.   Semi-connected patches

The escalating up case ($\alpha_{12} = 0$) is similar to the isolated case since patch two always dies out while patch 1 either reaches a coexistence equilibrium in the left plot or converges to a periodic solution in the right plot of Figure 3. Note that the parameters in the left (resp. right) plot for patch 1 of Figure 3 are identical to the left (resp. right) plot of Figure 2. Since we have shown that $R_{0, eu} = R_{0, 1}$, the global dynamics is not heavily impacted by the addition of uni-lateral host migration.

The cascading down model ($\alpha_{21} = 0$) presents three possible dynamics which depend on the value of $R_{0, cd} = \max(R_{0, 1}(1 - \alpha_{12}), R_{0, 2})$. The first case leads to extinction as it is shown in the left plot of Figure 4 where $R_{0, cd} = 0.75 < 1$. The second case leads to convergence to a coexistence equilibrium in the central plot of Figure 4 where $1 < R_{0, cd} = 3.41 < e^2$. The third case leads to convergence to a periodic solution as in the right plot of Figure 4 where $R_{0, cd} = 8.18 > e^2$.

7.3.   Interconnected patches

The interconnected case dynamics is similar to the escalating down model since in both patches ticks would either persist or die out. Though we have not shown full stability of this model, we observe similar behaviours as in the other cases. In particular, the threshold $R_{0, c}$ is key to determine the global stability of the model and can be computed using its definition. In the left plot of Figure 5, $R_{0, c} = 0.74 < 1$ guarantees tick extinction; in the central plot of Figure 5, dynamics lead to a coexistence equilibrium since $1 < R_{0, c} = 4.57 < e^2$ while in the right plot of Figure 5, tick dynamics converges to a periodic solution for $R_{0, c} = 8.02 > e^2$ and delay larger than a specific threshold. We also included a bifurcation diagram to show how the survival probability $\rho_1$ and the delay $\tau_1$ affect tick dynamics.

8.   Conclusions

We have introduced a spatial DDE model for tick demographics in a two-patch environment and shown how changes in environment favourability and tick movement on large mammals could affect the dynamics. Depending on the key parameter $R_{0, c}$ and on the delay parameters $\tau_1$ and $\tau_2$, we have shown three possible long-term behaviours of tick population including extinction, convergence to a coexistent solution and convergence to a periodic solution.

It is important to note that tick control measures can alter some of the model parameters used to differentiate three possible outcomes of tick population dynamics. In particular, patch-specific survival probabilities and development delays can be altered by habitat modification strategies including controlled burns. In addition, the mobility of large mammal hosts between patches can be modified by interventions such as deer fencing. So, our study here provides insights on how human interventions can change tick population dynamics. We recall that both thresholds are functions of the basic reproduction number in the isolated patches and of the mobility parameters. Using the closed form we obtained, we can see how these thresholds vary by changing these parameters. In particular, we conclude that

${\rm{(i)}}$ Increasing $R_{0, 1}$ helps coexistence since it decreases the threshold:

${\rm{(ii)}}$ Increasing $R_{0, 2}$ also helps coexistence, since

so if the first term of the product in the numerator is negative, coexistence is always possible (On the other hand, if this term is positive, then the threshold always decreases).

${\rm{(iii)}}$ If $\alpha_{21}$ is negligible with respect to $\alpha_{12}$ $(\alpha_{21} = o(1))$, then

Therefore $T_{coex}$ increases at first for small $\alpha_{12}$ and then undergoes a unique change of monotonicity at $\alpha_{12} = \sqrt{\frac{R_{0, 1}-1}{R_{0, 1}}}$. This implies that in order to facilitate coexistence, it is necessary to have a large movement of ticks from favorable to unfavorable environment. This threshold also depends on $R_{0, 1}$ which helps determine when the change of monotonicity occurs.

${\rm{(iv)}}$ If $\alpha_{12}$ is negligible with respect to $\alpha_{21}$ $(\alpha_{12} = o(1))$. Then increasing $\alpha_{21}$ does not facilitate coexistence since

There are a number of topics remaining for further studies. For example, the global stability of the coexistence equilibrium has not been resolved, though we suspect this is true and our numerical simulations also support this conjecture, in the case in which $e < R_{0, c} < e^2$. We suggest the idea that exponential ordering ^[33] coupled with the monotone dynamical systems theory can be used to establish the global convergence. We have not done much stability analysis for the bifurcated periodic solutions for (6.13), and the global continuation of these local Hopf branches also deserves additional attention ^[40]. Finally, in the real world setting, the environmental conditions are subject to seasonal temperature variations so a more realistic model capturing the tick population dynamics in the natural environment needs to be non-autonomous, at least periodic in time ^[41].

Acknowledgments

This research has been supported by the Natural Science Foundation of China (No. 12171074) (XZ), the Canada Research Chairs program and the Natural Sciences and Engineering Research Council of Canada (JW). We wanted to thank also Ning Yu for helping in the bifurcation analysis section.

Conflict of interest

The authors declare there is no conflict of interest.