Long-term transmission dynamics of tick-borne diseases involving seasonal variation and co-feeding transmission

Co-feeding is a mode of pathogen transmission for a wide range of tick-borne diseases where susceptible ticks can acquire infection from co-feeding with infected ticks on the same hosts. The significance of this transmission pathway is determined by the co-occurrence of ticks at different stages in the same season. Taking this into account, we formulate a system of differential equations with tick population dynamics and pathogen transmission dynamics highly regulated by the seasonal temperature variations. We examine the global dynamics of the model systems, and show that the two important ecological and epidemiological basic reproduction numbers can be used to fully characterize the long-term dynamics, and we link these two important threshold values to efficacy of co-feeding transmission.


Introduction
Incidences of tick-borne diseases (TBDs) such as Lyme disease or tick-borne encephalitis (TBE) have significantly increased in several parts of Europe and North America in the past decades. Climate conditions influence the transmission risk of TBDs since the life cycle of ticks depends on abiotic factors such as temperature, daylight length and humidity as well as biotic factors, for example the host abundance [1,2].
The TBD pathogens are delivered to a susceptible host during the blood-feeding of infected ticks. When this provokes a systemic infection in the host, susceptible ticks can subsequently get infected as they feed the infected host [26]. Recent discoveries suggest that pathogens may also be transmitted between co-feeding infected and uninfected ticks in the absence of the pathogen being established in the host [27,34]. This co-feeding transmission generally occurs between infected nymphs and susceptible larvae [26,34]. The significance of this mode of transmission depends on whether ticks in larval and nymphal stages are active in the same seasons of the year [26]. Co-feeding transmission is reported for many types of vector-borne pathogens and bridging hosts [6,27,34]. In particular, it is shown that Ixodes ricinis, a native European tick can efficiently transmit TBE virus and Lyme Borrelia bacterium via non-viraemic hosts [6,13].
The co-feeding transmission pathway is considered to enhance the transmission risk. Unlike systemic transmission, non-systemic transmission may occur immediately after the infectious tick taking blood meal and it dampens reduction of the transmission potential induced by natural mortality of hosts. In addition, more types of hosts involved in the disease transmission since some hosts which do not develop systemic infection can still serve as bridges for the co-feeding transmission [27]. Moreover, co-feeding transmission is observed in the hosts with immunity resulting in the increase of the number of transmission-competent hosts [14,26,27].
Several studies adapted a structured compartmental model as a framework to study the population dynamics of ticks and TBDs [5,18]. A basic reproduction number R 0 , a local stability threshold of such a dynamical system is often used to inform the risk factors [33]. Norman et al. [20] and Rosà and Pugliese [30] studied the effect of the host population structure and their densities on the tick population growth or TBD by identifying R 0 s of both the tick population and TBD transmission dynamics. Similarly, White at al. [36] studied a condition for the coexistence of two pathogens and Lou at al. [17] investigated how co-infection changes the fitness of pathogens.
In order to gain insights on which factors more significantly affect the long-term behaviours of tick population dynamics or the infection dynamics, some studied the stability or persistence in models for multistage populations [3,4,8,10,16,37,39]. Zhang and Wu [39] studied how the vector attachment and host grooming lead to bi-stability and nonlinear oscillations in the tick populations. Fan et al. [3] showed that a negative feedback from feeding to egg-laying adults leads to oscillations of tick population. A few studies qualitatively examined tick population or TBE transmission dynamics in a seasonal environment. For example, Lou et al. [16] and Wu and Wu [37] studied periodic systems describing population dynamics of ticks and TBD spread using the theory of monotone dynamical systems.
Some of the theoretical modelling studies included the co-feeding transmission in their models [7,21,31,36,38,40]. Using the pathogen basic reproduction number, Rosà et al. [31] showed that a co-feeding transmission enables the disease to persist even when hosts do not acquire and transmit the disease systemically. A similar result was obtained in [21]. In addition, some studies showed that co-feeding transmission increases the basic reproduction number or the prevalence in the study regions [19,28].
In this paper, we study the global dynamics of a proposed periodic system of ordinary differential equations describing the transmission of pathogens between ticks and hosts. The basic reproduction numbers for tick population growth and TBD spread are defined, and we will show that these two basic reproduction numbers determine the global longterm behaviours of solutions. We also perform some numerical simulations to see how co-feeding transmission changes the asymptotic behaviours of epidemiological system.

Tick population dynamics
We first introduce a model describing the dynamics of tick populations in various life stages: eggs (E), questing larvae (L q ), feeding larvae (L f ), questing nymphs (L q ), feeding nymphs (L f ), questing adults (A q ), feeding adults (A f ). Ticks in questing stages (L q , N q , A q ) move to feeding stages with the host-attaching rates α l (t), α n (t) and α a (t). Upon successful blood feeding, adult ticks will reproduce eggs with the rate b(t), while larval and nymphal ticks develop into the next stages with the development rates d l (t) and d n (t).
The mortality of ticks attached to hosts is density-dependent due to host immunity and grooming behaviours [15,25]. As in other modelling studies [17,22,37], we assume that ticks in each feeding stage experience density-dependent mortalities which is increasing with the number feeding ticks in the same stage. As an example, the mortality of feeding larvae is an increasing function of the average number of feeding larvae per host (L f /H). Questing ticks and eggs are assumed to undergo density-independent, constant mortality rates. Accounting for the seasonal dependence of the parameters, we describe the tick population dynamics with the following non-autonomous system of ordinary differential equations, where parameters are positive, bounded, p-periodic and the density-dependent mortalities μ fl (·), μ fn (·) and μ fa (·) are strictly increasing and continuously differentiable functions.
The tick population growth model with the same structure is studied in previous works and a proof of the following lemma can be found in paper [37].

Lemma 2.1:
The system (1) has a unique, non-negative and bounded solution for each initial value in R 7 + .

Disease transmission dynamics
We divide the tick and host population into those which are susceptible and infectious to TBD, with subscript 's' and 'i', respectively. All questing larvae are considered to be susceptible. In line with this assumption, we limit our focus on the tick-borne diseases with rare transovarial transmission, such as for Lyme disease or tick-borne encephalitis [23,24]. Upon the successful attachment to an infected host, the questing larvae and susceptible questing nymphs will be infected with probabilities β hl and β hn respectively, via the systemic infection route. In addition, we consider co-feeding infection, that is, when they feed the susceptible hosts or feed on infected hosts but escapes the systemic infection, the ticks may still get infected via non-systemic transmission with a probability δ(N fi /H). This probability is assumed to be increasing with respect to the average number of infected feeding nymphs per hosts. Finally, we assume that when the infected questing nymphs attach to a susceptible host, the host may be infected with a probability β nh . Therefore, we obtain a system describing the spread of TBD pathogens, where δ(x) = 1 − (1 − c) x . Here, c refers to the probability of non-systemic infection for susceptible and recovered ticks attached to a host attached with a single infected nymphs. We assume that all parameters are positive, bounded, p-periodic, 0 ≤ β hl ≤ 1, 0 ≤ β hn ≤ 1, 0 ≤ β nh ≤ 1 and 0 ≤ c < 1.
Let R n be n-dimensional vector space with Euclidean norm | · | n and C r : R → R n be the Banach space of p-periodic continuous functions. For

Lemma 2.2:
The system (2) has a unique, non-negative and bounded solution with initial values in R 12 + . Moreover, a feasible set Proof: With the assumption that the density-dependent moralities μ fl , μ fn and μ fa are continuously differentiable, it is easy to show that the vector field given by the right hand side of the system (2) is Lipschitz continuous on the bounded subsets of R 12 . Therefore, for a given initial value in R 12 , there exists a unique solution of the system (2) with the maximal interval of existence. It follows from Theorem 5.2.1 in [32] that solutions with initial values in R 12 + are non-negative. Let be a solution corresponding to a initial value in a feasible set X. Then, is a solution of the tick population dynamics (1) with initial value in R 7 + . By Lemma 2.1, y(t) is bounded. Since L fs + L fi is bounded and L fs , L fi are non-negative, L fs and L fi are bounded. Similarly, we can show that N q (t), We consider a feasible domain From Lemma 2.2, it follows that the feasible domain Z is positively invariant.

Tick population dynamics
Asymptotic behaviour of the tick population dynamics is studied in [37]. Herein, we introduce the results of this study. Basic reproduction number for the tick population, R tick is defined using the methods presented in [35]. The next generation operator L is defined as where F(t) and V(t) are periodic matrix valued functions representing new birth and transition between compartments taking forms as and Y(t, s), t ≥ s is the evolution operator of the systeṁ , the spectral radius of the next generation operator L.
which is globally stable with initial values in R 7 \ {0}.

The infection subsystem
When R tick > 1, there exists a unique positive p-periodic solution (5) which attracts every positive solution of the tick subsystem (4). Note that the system (3) possesses a disease-free p-periodic solution Consider the following limiting system arising at DFE(t), Linearizing (6) at an equilibrium (0, 0, 0, 0), we obtain the infected sub-system where F(t) and V(t) are time periodic matrix valued functions representing reproduction of new infections and transition between compartments taking forms as and The basic reproduction number is the spectral radius of the next generation operator defined with the above F(t) and V(t). Let R TBD := ρ(L), the spectral radius of L.
In the following, we prove global stability of the unique positive periodic solution following the method presented in [16]. (6). Then, P t is strongly monotone for all t ≥ 3p. i.e. for any t ≥ 3p and u, v ∈ R 4 , u < v implies P t (u) P t (v).

Lemma 3.2: Let P t be the solution map of the infected subsystem
Proof: We write the system (6) as We denote u(t, u 0 ) as the solution of (6) with initial value u 0 ∈ R 4 , u(t, w))). Then, X(0) = I and Note that It follows that for any given i, j ∈ {1, 2, 3, 4}, for all t ≥ 0. If there exists t 0 > 0 such that x ij (t 0 ) > 0, then it follows that x ij (t) > 0 for all t ≥ t 0 . Since x ii (0) = 1, we have x ii (t) > 0 for all t ≥ 0, i ∈ {1, 2, 3, 4}.

Characterizing the global dynamics of the full system
We can now characterize the global dynamics of the full system using the two basic reproductive numbers.

Theorem 3.4:
The global dynamics of the full system is completely determined by the three thresholds R tick and R TBD as follows: (1) If R tick ≤ 1, then zero is globally attractive for the system (3).
(3) If R tick > 1 and R TBD > 1, then there exists a positive periodic solution which is globally attractive for system (3) with respect to all positive solutions.
Proof: Let P be the Poincare map of the 11p-periodic system (3). In Lemma 2.2, we have shown that solutions of (3) are bounded and it follows that P is compact. Let ω = ω(x 0 ) be the omega limit set of P(x 0 ) for a given x 0 ∈ X. It follows from Lemma 2.1 of [9] that the omega limit set ω is an internally chain transitive set.
(2) Since R tick > 1, by Theorem 3.3, ω = {(l q (0),l f (0),n q (0),n f (0),â q (0),â f (0),m(0))} × ω 3 for some ω 3 ⊂ R 4 . Note that P(l q (0),l f (0),n q (0),n f (0),â q (0),â f (0),m(0), a) = (l q (0),l f (0),n q (0),n f (0),â q (0),â f (0),m(0), P 2 (a)) where P 2 is the Poincare map associated with the system (6). Since ω is an internally chain transitive set for P and it follows that ω 3 is an internally chain transitive set for P 2 . In the following, we prove that ω 3  When R TBD > 1, the continuity of the spectrum for matrices (Section 2.5.8 of [11]) implies that there exists δ > 0 such that the spectral radius of the Poincare map associated with the linearized system of the following system is greater than one: (10) Then, by the similar argument in the proof of Theorem 3.3, the system (10) has a positive periodic solution u * (t) which is globally attractive with respect to all positive solutions. Since R tick > 1, by Theorem 3.1, there exists some t 0 > 0 such that for all t > t 0 . By the comparison principle, which leads a contradiction to the assumption ω 3 = {0}. Since ω 3 = {0} and the positive periodic solution ee(t) is globally attractive for the nonzero solutions of the system (6), it follows that where W s (ee(0)) is the stable set for ee(0) with respect to the Poincare map P 2 . Then, by Theorem 3.1 of [9], By Theorem 3.1, the first 7 components of the globally attractive 11p-periodic solution is also p-periodic. Since the Poincare map P 2 of the 11p -periodic system has a unique positive fixed point, the fixed point is also a fixed point of the Poincare map associated to the p-periodic system. Therefore, the last three component of the globally attractive 11p-periodic solution is also p-periodic.

Discussions
In this study, we study the global dynamics of a system of differential equations describing transmission of pathogen among ticks and hosts. We identify two important ecological and epidemiological basic reproduction numbers and show that the two basic reproduction numbers can fully characterize the long-term dynamics. Co-feeding transmission formulated in our setting preserves the monotonicity of the dynamics of tick-host interaction for the pathogen transmission, and therefore the powerful theory of monotone dynamical systems theory can be applied to conclude the convergence to disease extinction or disease establishment status, characterized by periodic solutions due to the seasonal environmental condition variation. In our model, co-feeding transmission does not induce complex dynamics of the tick-host interaction. It is important to note, however, that co-feeding transmission can significantly alter the peak time/value and prevalence of disease transmissions. To illustrate this, we conduct a few numerical experiments. Figures 1 and 2 show the numerical solutions of the TBD transmission system with different values for cofeeding parameter, c. Table 1 show the parameters used in simulations. Most values are obtained from a modelling study of tick-borne encephalitis transmission [19], except for H/D (ratio of hosts for immature ticks to mature ticks) and T(t) (temperature at time t) which are assumed for the simulations. From the periodic attractors, we observe that the co-feeding parameter not only increases the peak of the periodic attractors of the infected populations, but it may also determine whether the disease persists or not as it is shown in other studies [7,27]. Regardless, contribution of co-feeding transmission in the perpetuation of a disease depends on the type of disease, hosts and tick abundance [12,29,34].  It should be mentioned that unlike our formulation of co-feeding transmission, feeding ticks are not equally distributed in hosts [27]. Some modelling studies on co-feeding transmission used a negative binomial distribution to describe the pattern of tick aggregation and studied the effect on disease dynamics of tick distribution [31,39]. Also, we assumed that the mortality of feeding ticks only depends on the tick population of the same stages but not feeding ticks in other stages. Considering that ticks in different life stages can cofeed on a same host, the density-dependent mortality of feeding ticks in some life-stage mortality rate of eggs 0.02618 μ ql , μ qn , μ qa mortality rates of questing 0.0068, 0.0034, 0.00136 larvae, nymphs, adults μ fl (l f ) mortality rate of feeding larvae 0.001428 + 0.01 · l f μ fn (n f ) mortality rate of feeding nymphs 0.000476 + 0.01 · n f μ fa (a f ) mortality rate of feeding adults 0.000408 + 0.1 · l f β hl , β hn , β nh transmission efficacy from hosts to larvae, 0.8, 0.8, 0.9 hosts to nymphs, nymphs to hosts c probability of non-systemic infection 0.4 when co-feeding with an infected nymph γ recovery rate of infected hosts 0. may depend on the populations of feeding ticks in other stages as well as its own population. In this case monotonicity of tick population dynamics is not guaranteed and it remains a topics for future research to examine to global dynamics of the corresponding model system.

Disclosure statement
No potential conflict of interest was reported by the author(s).