Mathematical model of Ehrlichia chaffeensis transmission dynamics in dogs

Abstract

Ehrlichia chaffeensis is a tick-borne disease transmitted by ticks to dogs. Few studies have mathematical modelled such tick-borne disease in dogs, and none have developed models that incorporate different ticks' developmental stages (discrete variable) as well as the duration of infection (continuous variable). In this study, we develop and analyze a model that considers these two structural variables using integrated semigroups theory. We address the well-posedness of the model and investigate the existence of steady states. The model exhibits a disease-free equilibrium and an endemic equilibrium. We calculate the reproduction number (${\mathcal{T}}_0$). We establish a necessary and sufficient condition for the bifurcation of an endemic equilibrium. Specifically, we demonstrate that a bifurcation, either backward or forward, can occur at ${\mathcal{T}}_0=1$, leading to the existence, or not, of an endemic equilibrium even when ${\mathcal{T}}_0<1$. Finally, numerical simulations are employed to illustrate these theoretical findings.

Keywords:

• Amblyomma americanum
• Ehrlichia chaffeensis
• age-structured model
• bifurcation analysis

MATHEMATICS SUBJECT CLASSIFICATIONs:

• 92-10
• 92D30

1. Introduction

Ehrlichiosis are tick-borne diseases caused by obligate intracellular rickettsia bacteria in the genera Ehrlichia [1–5]. These bacteria are classified within the group of the α-proteobacteria, order Rickettsiales, family Anaplasmataceae, genus Ehrlichia [3,6]. This genus consists of obligate intracellular Gram-negative bacteria [1,4,6] that mainly infect leukocytes (such as monocytes, macrophages, granulocytes and neutrophils), and endothelial cells in mammals, and salivary glands, intestinal epithelium, and hemolymph cells of ticks [1,7–12]. The genus Ehrlichia comprises of six recognized tick-transmitted species: E. canis, E. muris, E. chaffeensis, E. ewingii, E. minasensis, and E. ruminantium [3,11], and in recent times, other Ehrlichia species have been reported [7,8,11]. Going by our current knowledge, a large number of Ehrlichia species might not have been described [11].

The reservoir hosts for Ehrlichia species include numerous wild animals, as well as some domesticated species and livestock [6]. Some of these Ehrlichia species affect animals including pets such as cats and dogs [6,13], and a limited number have been know to infect humans [2,4,6,13–16]. Specifically, Ehrlichia canis, E. chaffeensis, E. ewingii, E. muris, and E. ruminantium are members of the genus Ehrlichia known to naturally infect various mammalian hosts such as cats, dogs, ruminants, and mice, and are responsible for emerging zoonoses in humans [7,8,10,17].

Ehrlichia are tick-borne diseases transmitted by ticks in the family Ixodidae. Rhipicephalus sanguineus (brown dog tick) and Ixodes ricinus are vectors for Ehrlichia canis [5,6]. E. canis can also be transmitted by Amblyomma cajennense, and experimentally by Dermacentor variabilis (American dog tick) [6]. In North America A. americanum (Lone Star tick) is the primary vector for both E. chaffeensis and E. ewingii [4,6,18,19]. ‘Outside North America, E. chaffeensis has been found in ticks in the genera Amblyomma, Haemaphysalis, Dermacentor, and Ixodes in Asia, in R. sanguineus in Cameroon, and in A. parvum in Argentina’ [6]. Haemaphysalis flava and Ixodes persulcatus complex ticks transmits E. muris [6]. E. ruminantium is the only known Ehrlichia species that infects cattle and it is found on the continent of Africa and a couple of Caribbean islands [8].

Ticks life-cycle span between 2-to-3 years depending on the species [20,21]. They go through four life stages namely egg, six-legged larva, eight-legged nymph, and adult. To survive each life stage after hatching from eggs, the ticks must take a blood meal from a host; but, most ticks will die if they are unable to find a host [20]. Ticks feed once on a host, then fall off and develop into the next stage. This feeding pattern creates pathways for diseases to be transmitted from hosts to hosts [22]. Ixodes scapularis (black-legged tick) life cycle generally lasts two years, while the life cycle of Amblyomma americanum (the lone star tick) is around three years long [21]. Some species like Rhipicephalus sanguineus (the brown dog tick) whose primary host is the domestic dog, prefer to feed on the same host during all its life stages [20,22]. Since R. sanguineus are endophilic, they are found inside houses and dog kennels [22].

Dogs are susceptible to infection with multiple Ehrlichia spp., including E. chaffeensis, E. ewingii, and E. canis [2,12,19,23,24]. In 2019, over 200,000 dogs tested positive for antibodies against Ehrlichia spp. within the United States out of 7,056,709 dogs tested, while over 1000 dogs tested positive out of 168,216 dogs tested in Canada [25]. According to Gettings et al. [25] the distribution of infected dogs follows the distribution of the related tick vectors. For instance, Amblyomma americanum is commonly found on dogs and people in the southeastern and southcentral United States [2]. In Beall et al. [2] the overall seroprevalence of E. canis, E. chaffeensis, and E. ewingii across the United States in 2012 was 0.8%, 2.8%, and 5.1%, respectively. The highest E. canis seroprevalence of 2.3% was found in Arkansas, Louisiana, Oklahoma, Tennessee, and Texas. E. chaffeensis seroreactivity was 6.6% in Arkansas, Kansas, Missouri, and Oklahoma (the central region), and 4.6% in Georgia, Maryland, North Carolina, South Carolina, Tennessee and Virginia (the southeast region). Seroreactivity of E. ewingii was highest in the central region with 14.6% value while the southeast region had a seroreactivity value of 5.9%.

Ehrlichia in dogs was discovered in the 1970s when military dogs were returning from the Vietnam war. They found this disease to be extremely severe in German Shepards, Doberman Pinschers, Belgium Malinois, and Siberian Huskies. Several studies including experimental researches and serological surveys have been carried out to understand Ehrlichiosis transmission in canine including Ehrlichia Chaffeensis [1,15,26–28]. We only found two quantitative studies modelling Ehrlichiosis in human [29,30]. Thus, in this study, we developed a mathematical model for Amblyomma americanum in the Great Plains and dogs infected with Ehrlichia Chaffeensis with the aim of understanding the quantitative properties of the model including the asymptotic dynamics of the disease in dogs across the Plains.

The rest of the paper is organized as follows: Section 2 gives the description of the model including the interactions between the ticks and dogs; Section 3 describes the results of the study. These results include the existence and uniqueness of solutions, the disease invasion process, and the bifurcation analysis. More precisely, we show that, depending on the sign of a constant ${\mathcal{C}}_{\mathrm{bif}}$ – referred to as the bifurcation parameter – which depends on model's parameters, a bifurcation occurs at ${\mathcal{T}}_0=1$ that is either a forward or backward bifurcation. In a forward bifurcation, which occurs when ${\mathcal{C}}_{\mathrm{bif}}<0$, there exists a unique endemic equilibrium if and only if ${\mathcal{T}}_0>1$. However, in a backward bifurcation, which occurs when ${\mathcal{C}}_{\mathrm{bif}}>0$, no endemic equilibrium exists for ${\mathcal{T}}_0\ll 1$ small enough; a unique endemic equilibrium exists if ${\mathcal{T}}_0>1$ while multiple equilibria exist when $0\ll {\mathcal{T}}_0<1$ close enough to 1. Finally in Section 4 we numerically illustrate such bifurcation results.

2. Model formulation

The transmission model of Ehrlichia Chaffeensis incorporates two subgroups: dogs and ticks. At any time t, the dog population is divided into susceptible $S_D(t)$, infected $i_D(t,a)$, and chronically infected $c_D(t,a)$, recovered $R_D(t)$. Here, the variable a represents the time since infection. Thus, the total number of dogs at time t is quantified by $\begin{array}{r}{N}_D(t)=S_D(t)+{\int }_0^{\mathrm{\infty }}{i}_D(t,a)\mathrm{d}a+{\int }_0^{\mathrm{\infty }}{c}_D(t,a)\mathrm{d}a+R_D(t).\end{array}$The ticks population is structured into several stages: eggs (E), larva (L), nymph (N), and adult (A). We denote by $\mathcal{S}=\{E,L,N,A\}$ the set of ticks stages. Let $S_{Tk}(t)$ be the number of susceptible ticks of stage $K\in \mathcal{S}$ at time t. We also denote by $i_{Tk}(t,a)$ the number of infected ticks of stage $K\in \mathcal{S}\setminus \{E\}$ at time t and which are infected since time a. The model proposed assumes that there are no infected eggs. Therefore the total number of ticks of stage $k\in \mathcal{S}\setminus \{E\}$ at time t is given by (1) $\begin{array}{r}{N}_{Tk}(t)=S_{Tk}(t)+{\int }_0^{\mathrm{\infty }}{i}_{Tk}(t,a)\mathrm{d}a.\end{array}$(1) The dogs-ticks infection life cycle is shown in Figure 1.

Figure 1. Flow chart of the model.

Dogs' population dynamic. At any time t, infected ticks (which are infected since time a) induce an infection within the dogs' population through the force of infection ${\lambda }_T(t,a)$, such that (2) $\begin{array}{r}{\lambda }_T(t,a)={\beta }_T(a)\sum _{K\in \mathcal{S}\setminus \{E\}}{i}_{Tk}(t,a),\end{array}$(2) where ${\beta }_T$ denotes the infectivity of infected ticks. Therefore, newly infected dogs are given by $\frac{S_D(t)}{N_D(t)}{\int }_0^{\mathrm{\infty }}{\lambda }_T(t,a)\mathrm{d}a$ and $ϵ\frac{R_D(t)}{N_D(t)}{\int }_0^{\mathrm{\infty }}{\lambda }_T(t,a)\mathrm{d}a$, where ϵ is constant parameter accounting for the relative disease transmission to previously recovered dogs. Note that ϵ can be considered age-dependent ($ϵ=ϵ(a)$) with no more difficulties in the analysis proposed here. Thus, the dogs population dynamic is then described by the below system (3) $\begin{array}{r}\{\begin{array}{l}{\displaystyle \frac{\mathrm{d}{S}_D(t)}{\mathrm{d}t}={\theta }_D-\frac{S_D(t)}{N_D(t)}{\int }_0^{\mathrm{\infty }}{\lambda }_T(t,a)\mathrm{d}a-{\mu }_D{S}_D(t),}\\ i_D(t,0)={\displaystyle \frac{S_D(t)+ϵR_D(t)}{N_D(t)}{\int }_0^{\mathrm{\infty }}{\lambda }_T(t,a)\mathrm{d}a,}\\ {\displaystyle \frac{\mathrm{\partial }{i}_D(t,a)}{\mathrm{\partial }t}+\frac{\mathrm{\partial }{i}_D(t,a)}{\mathrm{\partial }a}=-({\mu }_D+\gamma (a)+{\delta }_D(a)+{\nu }_D(a))i_D(t,a),}\\ {\displaystyle \frac{\mathrm{\partial }{c}_D(t,a)}{\mathrm{\partial }t}+\frac{\mathrm{\partial }{c}_D(t,a)}{\mathrm{\partial }a}=-({\mu }_D+{\delta }_C(a))c_D(t,a)}\\ c_D(t,0)={\displaystyle {\int }_0^{\mathrm{\infty }}{\nu }_D(a)i_D(t,a)\mathrm{d}a,}\\ {\displaystyle \frac{\mathrm{d}{R}_D(t)}{\mathrm{d}t}={\int }_0^{\mathrm{\infty }}\gamma (a)i_D(t,a)\mathrm{d}a-ϵ\frac{R_D(t)}{N_D(t)}{\int }_0^{\mathrm{\infty }}{\lambda }_T(t,a)\mathrm{d}a-{\mu }_D{R}_D(t).}\end{array}\end{array}$(3) In the above model, susceptible dogs are recruited at rate ${\theta }_D$ and all dogs die naturally at rate ${\mu }_D$. Infected and chronically infected dogs have a disease-induced mortality rate ${\delta }_D$ and ${\delta }_C$. Infected dogs progress to a chronic infection at rate ${\nu }_D$. Finally, only non-chronic infections are assumed to recover from the infection at rate γ.

Ticks' population dynamic. At time t, infected dogs (chronic or not) induce an infection within the ticks population through the force of infection ${\lambda }_D(t)$, such that $\begin{array}{r}{\lambda }_D(t)={\int }_0^{\mathrm{\infty }}[{\beta }_D^1(a)i_D(t,a)+{\beta }_D^2(a)c_D(t,a)]\mathrm{d}a,\end{array}$where ${\beta }_D^j(a)$ denotes the infectivity of infected and chronically infected dogs a-time post infection. Therefore, newly infected ticks of stage $k\in \mathcal{S}\setminus \{E\}$ are given by $\frac{S_{Tk}(t)}{N_D(t)}{\lambda }_D(t)$. The ticks population dynamic is then described by the system below (4) $\begin{array}{r}\{\begin{array}{l}{\displaystyle \frac{\mathrm{d}{S}_{TE}(t)}{\mathrm{d}t}=r_E(1-\frac{S_{TE}(t)}{K})N_{TA}(t)-({\alpha }_E+{\mu }_E)S_{TE}(t),}\\ {\displaystyle \frac{\mathrm{d}{S}_{TL}(t)}{\mathrm{d}t}={\alpha }_E{S}_{TE}(t)-\frac{S_{TL}(t)}{N_D(t)}{\lambda }_D(t)-({\alpha }_L+{\mu }_L)S_{TL}(t),}\\ {\displaystyle \frac{\mathrm{d}{S}_{TN}(t)}{\mathrm{d}t}={\alpha }_L{S}_{TL}(t)-\frac{S_{TN}(t)}{N_D(t)}{\lambda }_D(t)-({\alpha }_N+{\mu }_N)S_{TN}(t),}\\ {\displaystyle \frac{\mathrm{d}{S}_{TA}(t)}{\mathrm{d}t}={\alpha }_N{S}_{TN}(t)-\frac{S_{TA}(t)}{N_D(t)}{\lambda }_D(t)-{\mu }_A{S}_{TA}(t),}\\ \text{and for }k\in \mathcal{S}\setminus \{E\},\\ i_{Tk}(t,0)={\displaystyle \frac{S_{Tk}(t)}{N_D(t)}{\lambda }_D(t),}\\ {\displaystyle \frac{\mathrm{\partial }{i}_{Tk}(t,a)}{\mathrm{\partial }t}+\frac{\mathrm{\partial }{i}_{Tk}(t,a)}{\mathrm{\partial }a}=-({\mu }_k+{\alpha }_k)i_{Tk}(t,a)}\end{array}\end{array}$(4) with ${\alpha }_A=0$. In the above system of the ticks dynamic, the eggs' production rate at time t is $r_E(1-\frac{S_{TE}(t)}{K})N_{TA}(t)$, where $r_E$ is the number of eggs produced by adult ticks $N_{TA}$, and the parameter K is the eggs' carrying capacity. Parameters ${\alpha }_k$,s are ticks progression rates from the eggs' to adult' stage. The death rate of ticks at each k-stage is ${\mu }_k$. The notations of all variables and parameters are summarized in Table 1.

Table 1. State variables and parameters used for simulations.

Category	Description	Value (Unit) [Ref]
Variables
t	Time	day
a	Time since infection	day
State Variables
$S_D(t)$	Susceptible Dogs
$i_D(t,a)$	Infected Dogs
$c_D(t,a)$	Chronically infected Dogs
$R_D(t)$	Recovered Dogs
$S_{TE}(t)$	Number of Eggs
$S_{Tk}(t)$	Susceptible ticks at stage k
$k\in \{L,N,A\}$	L for Larva, N for Nymph, and A for Adult
$i_{Tk}(t,a)$	Infected ticks at stage k
$k\in \{L,N,A\}$	L for Larva, N for Nymph, and A for Adult
Parameters
${\theta }_D$	Birth rate of dogs	$\mathrm{70,000}$ (day${}^{-1}$) [31]
${\mu }_D$	Death rate of dogs	0.00027397 (day${}^{-1}$) [32]
ϵ	Rate of reinfection in dogs	0.00444444 (day${}^{-1}$) [33]
$\gamma (a)$	Rate of infected to recovered dogs	0.04761905 (day${}^{-1}$) [34]
${\nu }_D(a)$	Rate of acute to chronic infection in dogs	0.04761905 (day${}^{-1}$) [34]
${\delta }_D(a)$	Acute infection death rate in dogs	Variable (day${}^{-1}$)
${\delta }_C(a)$	Chronic infection death rate in dogs	Variable (day${}^{-1}$)
${\beta }_D^1(a)$	Disease transmission probability from	Variable (No Unit)
	infected dogs to ticks
${\beta }_D^2(a)$	Disease transmission probability from	Variable (No Unit)
	chronically infected dogs to ticks
${\beta }_T(a)$	Disease transmission probability from ticks to dogs	Variable (No Unit)
$r_E$	Ticks egg laying rate	Variable (day${}^{-1}$)
K	Ticks carrying capacity	$\mathrm{1,000,000}$ (day${}^{-1}$) [Assumed]
${\alpha }_E$	Maturation rate from tick eggs to larvae	0.0243902 (day${}^{-1}$) [35]
${\alpha }_L$	Maturation rate from tick larvae to nymphs	0.00273973 (day${}^{-1}$) [35]
${\alpha }_N$	Maturation rate from tick nymphs to adult	0.0037037 (day${}^{-1}$) [35]
${\mu }_E$	Death rate of eggs	0.008 (day${}^{-1}$) [36]
${\mu }_k$	Death rate of ticks at stage k	0.003 (day${}^{-1}$) [Assumed]
$k\in \{L,N,A\}$	L for Larva, N for Nymph, and A for Adult

3. Main results

This section is devoted to the main results of this paper. Overall, our results are based on the following assumption on the parameters of system (3(3) $\begin{array}{r}\{\begin{array}{l}{\displaystyle \frac{\mathrm{d}{S}_D(t)}{\mathrm{d}t}={\theta }_D-\frac{S_D(t)}{N_D(t)}{\int }_0^{\mathrm{\infty }}{\lambda }_T(t,a)\mathrm{d}a-{\mu }_D{S}_D(t),}\\ i_D(t,0)={\displaystyle \frac{S_D(t)+ϵR_D(t)}{N_D(t)}{\int }_0^{\mathrm{\infty }}{\lambda }_T(t,a)\mathrm{d}a,}\\ {\displaystyle \frac{\mathrm{\partial }{i}_D(t,a)}{\mathrm{\partial }t}+\frac{\mathrm{\partial }{i}_D(t,a)}{\mathrm{\partial }a}=-({\mu }_D+\gamma (a)+{\delta }_D(a)+{\nu }_D(a))i_D(t,a),}\\ {\displaystyle \frac{\mathrm{\partial }{c}_D(t,a)}{\mathrm{\partial }t}+\frac{\mathrm{\partial }{c}_D(t,a)}{\mathrm{\partial }a}=-({\mu }_D+{\delta }_C(a))c_D(t,a)}\\ c_D(t,0)={\displaystyle {\int }_0^{\mathrm{\infty }}{\nu }_D(a)i_D(t,a)\mathrm{d}a,}\\ {\displaystyle \frac{\mathrm{d}{R}_D(t)}{\mathrm{d}t}={\int }_0^{\mathrm{\infty }}\gamma (a)i_D(t,a)\mathrm{d}a-ϵ\frac{R_D(t)}{N_D(t)}{\int }_0^{\mathrm{\infty }}{\lambda }_T(t,a)\mathrm{d}a-{\mu }_D{R}_D(t).}\end{array}\end{array}$(3) )–(4(4) $\begin{array}{r}\{\begin{array}{l}{\displaystyle \frac{\mathrm{d}{S}_{TE}(t)}{\mathrm{d}t}=r_E(1-\frac{S_{TE}(t)}{K})N_{TA}(t)-({\alpha }_E+{\mu }_E)S_{TE}(t),}\\ {\displaystyle \frac{\mathrm{d}{S}_{TL}(t)}{\mathrm{d}t}={\alpha }_E{S}_{TE}(t)-\frac{S_{TL}(t)}{N_D(t)}{\lambda }_D(t)-({\alpha }_L+{\mu }_L)S_{TL}(t),}\\ {\displaystyle \frac{\mathrm{d}{S}_{TN}(t)}{\mathrm{d}t}={\alpha }_L{S}_{TL}(t)-\frac{S_{TN}(t)}{N_D(t)}{\lambda }_D(t)-({\alpha }_N+{\mu }_N)S_{TN}(t),}\\ {\displaystyle \frac{\mathrm{d}{S}_{TA}(t)}{\mathrm{d}t}={\alpha }_N{S}_{TN}(t)-\frac{S_{TA}(t)}{N_D(t)}{\lambda }_D(t)-{\mu }_A{S}_{TA}(t),}\\ \text{and for }k\in \mathcal{S}\setminus \{E\},\\ i_{Tk}(t,0)={\displaystyle \frac{S_{Tk}(t)}{N_D(t)}{\lambda }_D(t),}\\ {\displaystyle \frac{\mathrm{\partial }{i}_{Tk}(t,a)}{\mathrm{\partial }t}+\frac{\mathrm{\partial }{i}_{Tk}(t,a)}{\mathrm{\partial }a}=-({\mu }_k+{\alpha }_k)i_{Tk}(t,a)}\end{array}\end{array}$(4) ). More precisely, we assume that

Assumption 3.1

(1)	The recruitment rate ${\theta }_D$ of susceptible dogs and the natural death rate ${\mu }_D$ are positive constants. The ticks' stage transition rates ${\alpha }_k$,s and death rates ${\mu }_k$,s are positive constants. The parameter ϵ, $r_E$ and the carrying capacity K are positive constants;
(2)	Parameters $\gamma ,{\delta }_D,{\nu }_D,{\delta }_C$ and the transmission rates ${\beta }_T,{\beta }_D^k$, k = 1, 2, belong in $L_{+}^{\mathrm{\infty }}((0,+\mathrm{\infty }),\mathbb{R})\setminus \{0_{L^{\mathrm{\infty }}}\}$;
(3)	The initial condition is such that $S_D(0)=S_{D0}>0$, $R_D(0)=R_{D0}\ge 0$, $i_D(0,\cdot )=i_{D0}\in L_{+}^1((0,+\mathrm{\infty }),\mathbb{R})$, $c_D(0,\cdot )=c_{D0}\in L_{+}^1((0,+\mathrm{\infty }),\mathbb{R})$, $S_{Tk0}(0)=S_{Tk0}\ge 0$, and $i_{Tk}(0,\cdot )=i_{Tk0}\in L_{+}^1((0,+\mathrm{\infty }),\mathbb{R})$ for $k\in \mathcal{S}$.

3.1. Existence and uniqueness of nonnegative solution

In this section, we state the results concerning the existence of a globally defined nonnegative solution to system (3(3) $\begin{array}{r}\{\begin{array}{l}{\displaystyle \frac{\mathrm{d}{S}_D(t)}{\mathrm{d}t}={\theta }_D-\frac{S_D(t)}{N_D(t)}{\int }_0^{\mathrm{\infty }}{\lambda }_T(t,a)\mathrm{d}a-{\mu }_D{S}_D(t),}\\ i_D(t,0)={\displaystyle \frac{S_D(t)+ϵR_D(t)}{N_D(t)}{\int }_0^{\mathrm{\infty }}{\lambda }_T(t,a)\mathrm{d}a,}\\ {\displaystyle \frac{\mathrm{\partial }{i}_D(t,a)}{\mathrm{\partial }t}+\frac{\mathrm{\partial }{i}_D(t,a)}{\mathrm{\partial }a}=-({\mu }_D+\gamma (a)+{\delta }_D(a)+{\nu }_D(a))i_D(t,a),}\\ {\displaystyle \frac{\mathrm{\partial }{c}_D(t,a)}{\mathrm{\partial }t}+\frac{\mathrm{\partial }{c}_D(t,a)}{\mathrm{\partial }a}=-({\mu }_D+{\delta }_C(a))c_D(t,a)}\\ c_D(t,0)={\displaystyle {\int }_0^{\mathrm{\infty }}{\nu }_D(a)i_D(t,a)\mathrm{d}a,}\\ {\displaystyle \frac{\mathrm{d}{R}_D(t)}{\mathrm{d}t}={\int }_0^{\mathrm{\infty }}\gamma (a)i_D(t,a)\mathrm{d}a-ϵ\frac{R_D(t)}{N_D(t)}{\int }_0^{\mathrm{\infty }}{\lambda }_T(t,a)\mathrm{d}a-{\mu }_D{R}_D(t).}\end{array}\end{array}$(3) )–(4(4) $\begin{array}{r}\{\begin{array}{l}{\displaystyle \frac{\mathrm{d}{S}_{TE}(t)}{\mathrm{d}t}=r_E(1-\frac{S_{TE}(t)}{K})N_{TA}(t)-({\alpha }_E+{\mu }_E)S_{TE}(t),}\\ {\displaystyle \frac{\mathrm{d}{S}_{TL}(t)}{\mathrm{d}t}={\alpha }_E{S}_{TE}(t)-\frac{S_{TL}(t)}{N_D(t)}{\lambda }_D(t)-({\alpha }_L+{\mu }_L)S_{TL}(t),}\\ {\displaystyle \frac{\mathrm{d}{S}_{TN}(t)}{\mathrm{d}t}={\alpha }_L{S}_{TL}(t)-\frac{S_{TN}(t)}{N_D(t)}{\lambda }_D(t)-({\alpha }_N+{\mu }_N)S_{TN}(t),}\\ {\displaystyle \frac{\mathrm{d}{S}_{TA}(t)}{\mathrm{d}t}={\alpha }_N{S}_{TN}(t)-\frac{S_{TA}(t)}{N_D(t)}{\lambda }_D(t)-{\mu }_A{S}_{TA}(t),}\\ \text{and for }k\in \mathcal{S}\setminus \{E\},\\ i_{Tk}(t,0)={\displaystyle \frac{S_{Tk}(t)}{N_D(t)}{\lambda }_D(t),}\\ {\displaystyle \frac{\mathrm{\partial }{i}_{Tk}(t,a)}{\mathrm{\partial }t}+\frac{\mathrm{\partial }{i}_{Tk}(t,a)}{\mathrm{\partial }a}=-({\mu }_k+{\alpha }_k)i_{Tk}(t,a)}\end{array}\end{array}$(4) ).

In order to state our main results, we will rewrite the system in a more appropriate equivalent form. To do this, denote by ${\mathbf{e}}_k,k\in \{1,2,3\}$ and ${\mathbf{v}}_k$, $k\in \{1,2\}$ the canonical basis of ${\mathbb{R}}^3$ and ${\mathbb{R}}^2$, respectively. Thus, we can consider the states variables in a vector form by setting for every $t\ge 0$ and $a\ge 0$ $\begin{array}{r}\{\begin{array}{l}{\mathit{S}}_T(t)=S_{TL}(t){\mathbf{e}}_1+S_{TN}(t){\mathbf{e}}_2+S_{TA}(t){\mathbf{e}}_3,\\ {\mathit{i}}_T(t,a)=i_{TL}(t,a){\mathbf{e}}_1+i_{TN}(t,a){\mathbf{e}}_2+i_{TA}(t,a){\mathbf{e}}_3,\\ {\mathit{i}}_D(t,a)=i_D(t,a){\mathbf{v}}_1+c_D(t,a){\mathbf{v}}_2.\end{array}\end{array}$Moreover, we define the vector ${\mathit{\beta }}_D(a)$ of transmission rates with components ${\beta }_D^1(a)$ and ${\beta }_D^2(a)$, i.e. $\begin{array}{r}{\mathit{\beta }}_D(a)={\beta }_D^1(a){\mathbf{v}}_1+{\beta }_D^2(a){\mathbf{v}}_2,\end{array}$and define the vector $\mathbf{1}={\mathbf{e}}_1+{\mathbf{e}}_2+{\mathbf{e}}_3$. Let $\mathbf{M}$ be ticks stage progression matrix from the larval stage to the adult stage, i.e. (5) $\begin{array}{r}\mathbf{M}:=\left(\begin{array}{ccc}{\mu }_L+{\alpha }_L& 0& 0\\ -{\alpha }_L& {\mu }_N+{\alpha }_N& 0\\ 0& -{\alpha }_N& {\mu }_A\end{array}\right).\end{array}$(5) Therefore, using the notation $⟨\cdot ,\cdot ⟩$ to denotes the inner product in ${\mathbb{R}}^2$ and ${\mathbb{R}}^3$ and setting for all $a\ge 0$ (6) $\begin{array}{r}\mathit{\theta }(a):=\left(\begin{array}{cc}\gamma (a)+{\delta }_D(a)+{\nu }_D(a)& 0\\ 0& {\delta }_C(a)\end{array}\right),\end{array}$(6) Model (3(3) $\begin{array}{r}\{\begin{array}{l}{\displaystyle \frac{\mathrm{d}{S}_D(t)}{\mathrm{d}t}={\theta }_D-\frac{S_D(t)}{N_D(t)}{\int }_0^{\mathrm{\infty }}{\lambda }_T(t,a)\mathrm{d}a-{\mu }_D{S}_D(t),}\\ i_D(t,0)={\displaystyle \frac{S_D(t)+ϵR_D(t)}{N_D(t)}{\int }_0^{\mathrm{\infty }}{\lambda }_T(t,a)\mathrm{d}a,}\\ {\displaystyle \frac{\mathrm{\partial }{i}_D(t,a)}{\mathrm{\partial }t}+\frac{\mathrm{\partial }{i}_D(t,a)}{\mathrm{\partial }a}=-({\mu }_D+\gamma (a)+{\delta }_D(a)+{\nu }_D(a))i_D(t,a),}\\ {\displaystyle \frac{\mathrm{\partial }{c}_D(t,a)}{\mathrm{\partial }t}+\frac{\mathrm{\partial }{c}_D(t,a)}{\mathrm{\partial }a}=-({\mu }_D+{\delta }_C(a))c_D(t,a)}\\ c_D(t,0)={\displaystyle {\int }_0^{\mathrm{\infty }}{\nu }_D(a)i_D(t,a)\mathrm{d}a,}\\ {\displaystyle \frac{\mathrm{d}{R}_D(t)}{\mathrm{d}t}={\int }_0^{\mathrm{\infty }}\gamma (a)i_D(t,a)\mathrm{d}a-ϵ\frac{R_D(t)}{N_D(t)}{\int }_0^{\mathrm{\infty }}{\lambda }_T(t,a)\mathrm{d}a-{\mu }_D{R}_D(t).}\end{array}\end{array}$(3) )–(4(4) $\begin{array}{r}\{\begin{array}{l}{\displaystyle \frac{\mathrm{d}{S}_{TE}(t)}{\mathrm{d}t}=r_E(1-\frac{S_{TE}(t)}{K})N_{TA}(t)-({\alpha }_E+{\mu }_E)S_{TE}(t),}\\ {\displaystyle \frac{\mathrm{d}{S}_{TL}(t)}{\mathrm{d}t}={\alpha }_E{S}_{TE}(t)-\frac{S_{TL}(t)}{N_D(t)}{\lambda }_D(t)-({\alpha }_L+{\mu }_L)S_{TL}(t),}\\ {\displaystyle \frac{\mathrm{d}{S}_{TN}(t)}{\mathrm{d}t}={\alpha }_L{S}_{TL}(t)-\frac{S_{TN}(t)}{N_D(t)}{\lambda }_D(t)-({\alpha }_N+{\mu }_N)S_{TN}(t),}\\ {\displaystyle \frac{\mathrm{d}{S}_{TA}(t)}{\mathrm{d}t}={\alpha }_N{S}_{TN}(t)-\frac{S_{TA}(t)}{N_D(t)}{\lambda }_D(t)-{\mu }_A{S}_{TA}(t),}\\ \text{and for }k\in \mathcal{S}\setminus \{E\},\\ i_{Tk}(t,0)={\displaystyle \frac{S_{Tk}(t)}{N_D(t)}{\lambda }_D(t),}\\ {\displaystyle \frac{\mathrm{\partial }{i}_{Tk}(t,a)}{\mathrm{\partial }t}+\frac{\mathrm{\partial }{i}_{Tk}(t,a)}{\mathrm{\partial }a}=-({\mu }_k+{\alpha }_k)i_{Tk}(t,a)}\end{array}\end{array}$(4) ) takes the following form (7) $\begin{array}{r}\{\begin{array}{l}{\displaystyle \frac{\mathrm{d}{S}_D(t)}{\mathrm{d}t}={\theta }_D-\frac{S_D(t)}{N_D(t)}{\int }_0^{\mathrm{\infty }}{\beta }_T(a)⟨\mathbf{1},{\mathit{i}}_T(t,a)⟩\mathrm{d}a-{\mu }_D{S}_D(t),}\\ {\displaystyle \frac{\mathrm{\partial }{\mathit{i}}_D(t,a)}{\mathrm{\partial }t}+\frac{\mathrm{\partial }{\mathit{i}}_D(t,a)}{\mathrm{\partial }a}=-({\mu }_D+\mathit{\theta }(a)){\mathit{i}}_D(t,a),}\\ {\mathit{i}}_D(t,0)={\displaystyle \frac{S_D(t)+ϵR_D(t)}{N_D(t)}{\int }_0^{\mathrm{\infty }}{\beta }_T(a)⟨\mathbf{1},{\mathit{i}}_T(t,a)⟩\mathrm{d}a{\mathbf{v}}_1}\\ {\displaystyle +{\int }_0^{\mathrm{\infty }}{\nu }_D(a)⟨{\mathbf{v}}_1,{\mathit{i}}_D(t,a)⟩\mathrm{d}a{\mathbf{v}}_2,}\\ {\displaystyle \frac{\mathrm{d}{R}_D(t)}{\mathrm{d}t}={\int }_0^{\mathrm{\infty }}\gamma (a)⟨{\mathbf{v}}_1,{\mathit{i}}_D(t,a)⟩\mathrm{d}a-ϵ\frac{R_D(t)}{N_D(t)}{\int }_0^{\mathrm{\infty }}{\beta }_T(a)⟨\mathbf{1},{\mathit{i}}_T(t,a)⟩\mathrm{d}a-{\mu }_D{R}_D(t),}\end{array}\end{array}$(7) and (8) $\begin{array}{r}\{\begin{array}{l}{\displaystyle \frac{\mathrm{d}{S}_{TE}(t)}{\mathrm{d}t}=r_E(1-\frac{S_{TE}(t)}{K}){\displaystyle (⟨{\mathit{S}}_T(t),{\mathbf{e}}_3⟩+{\int }_0^{\mathrm{\infty }}⟨{\mathit{i}}_T(t,a),{\mathbf{e}}_3⟩\mathrm{d}a)-({\alpha }_E+{\mu }_E)S_{TE}(t),}}\\ {\displaystyle \frac{\mathrm{d}{\mathit{S}}_T(t)}{\mathrm{d}t}={\alpha }_E{S}_{TE}(t){\mathbf{e}}_1-\mathbf{M}{\mathit{S}}_T(t)-\frac{{\mathit{S}}_T(t)}{N_D(t)}{\int }_0^{\mathrm{\infty }}⟨{\mathit{\beta }}_D(a),{\mathit{i}}_D(t,a)⟩\mathrm{d}a,}\\ {\displaystyle \frac{\mathrm{\partial }{\mathit{i}}_T(t,a)}{\mathrm{\partial }t}+\frac{\mathrm{\partial }{\mathit{i}}_T(t,a)}{\mathrm{\partial }a}=-\mathbf{M}{\mathit{i}}_T(t,a),}\\ {\mathit{i}}_T(t,0)={\displaystyle \frac{{\mathit{S}}_T(t)}{N_D(t)}{\int }_0^{\mathrm{\infty }}⟨{\mathit{\beta }}_D(a),{\mathit{i}}_D(t,a)⟩\mathrm{d}a,}\end{array}\end{array}$(8) subject to the initial condition $S_D(0)=S_{D0}>0$, ${\mathit{i}}_D(0,\cdot )={\mathit{i}}_{D0}\in L_{+}^1((0,+\mathrm{\infty }),{\mathbb{R}}^2)$, $R_D(0)=R_{D0}\in {\mathbb{R}}_{+}$, $S_{TE}(0)=S_{TE0}\in {\mathbb{R}}_{+}$, ${\mathit{i}}_D(0,\cdot )={\mathit{i}}_{D0}\in L_{+}^1((0,+\mathrm{\infty }),{\mathbb{R}}^2)$, ${\mathit{S}}_T(0)={\mathit{S}}_{T0}\in {\mathbb{R}}_{+}^3$, and ${\mathit{i}}_T(0,\cdot )={\mathit{i}}_{T0}\in L_{+}^1((0,+\mathrm{\infty }),{\mathbb{R}}^3)$. The following result then concerns the existence and uniqueness of nonnegative solutions to (7(7) $\begin{array}{r}\{\begin{array}{l}{\displaystyle \frac{\mathrm{d}{S}_D(t)}{\mathrm{d}t}={\theta }_D-\frac{S_D(t)}{N_D(t)}{\int }_0^{\mathrm{\infty }}{\beta }_T(a)⟨\mathbf{1},{\mathit{i}}_T(t,a)⟩\mathrm{d}a-{\mu }_D{S}_D(t),}\\ {\displaystyle \frac{\mathrm{\partial }{\mathit{i}}_D(t,a)}{\mathrm{\partial }t}+\frac{\mathrm{\partial }{\mathit{i}}_D(t,a)}{\mathrm{\partial }a}=-({\mu }_D+\mathit{\theta }(a)){\mathit{i}}_D(t,a),}\\ {\mathit{i}}_D(t,0)={\displaystyle \frac{S_D(t)+ϵR_D(t)}{N_D(t)}{\int }_0^{\mathrm{\infty }}{\beta }_T(a)⟨\mathbf{1},{\mathit{i}}_T(t,a)⟩\mathrm{d}a{\mathbf{v}}_1}\\ {\displaystyle +{\int }_0^{\mathrm{\infty }}{\nu }_D(a)⟨{\mathbf{v}}_1,{\mathit{i}}_D(t,a)⟩\mathrm{d}a{\mathbf{v}}_2,}\\ {\displaystyle \frac{\mathrm{d}{R}_D(t)}{\mathrm{d}t}={\int }_0^{\mathrm{\infty }}\gamma (a)⟨{\mathbf{v}}_1,{\mathit{i}}_D(t,a)⟩\mathrm{d}a-ϵ\frac{R_D(t)}{N_D(t)}{\int }_0^{\mathrm{\infty }}{\beta }_T(a)⟨\mathbf{1},{\mathit{i}}_T(t,a)⟩\mathrm{d}a-{\mu }_D{R}_D(t),}\end{array}\end{array}$(7) )–(8(8) $\begin{array}{r}\{\begin{array}{l}{\displaystyle \frac{\mathrm{d}{S}_{TE}(t)}{\mathrm{d}t}=r_E(1-\frac{S_{TE}(t)}{K}){\displaystyle (⟨{\mathit{S}}_T(t),{\mathbf{e}}_3⟩+{\int }_0^{\mathrm{\infty }}⟨{\mathit{i}}_T(t,a),{\mathbf{e}}_3⟩\mathrm{d}a)-({\alpha }_E+{\mu }_E)S_{TE}(t),}}\\ {\displaystyle \frac{\mathrm{d}{\mathit{S}}_T(t)}{\mathrm{d}t}={\alpha }_E{S}_{TE}(t){\mathbf{e}}_1-\mathbf{M}{\mathit{S}}_T(t)-\frac{{\mathit{S}}_T(t)}{N_D(t)}{\int }_0^{\mathrm{\infty }}⟨{\mathit{\beta }}_D(a),{\mathit{i}}_D(t,a)⟩\mathrm{d}a,}\\ {\displaystyle \frac{\mathrm{\partial }{\mathit{i}}_T(t,a)}{\mathrm{\partial }t}+\frac{\mathrm{\partial }{\mathit{i}}_T(t,a)}{\mathrm{\partial }a}=-\mathbf{M}{\mathit{i}}_T(t,a),}\\ {\mathit{i}}_T(t,0)={\displaystyle \frac{{\mathit{S}}_T(t)}{N_D(t)}{\int }_0^{\mathrm{\infty }}⟨{\mathit{\beta }}_D(a),{\mathit{i}}_D(t,a)⟩\mathrm{d}a,}\end{array}\end{array}$(8) ).

Theorem 3.2

Let Assumption 3.1 be satisfied. Then there exists a unique continuous globally defined integrated solution to (7(7) $\begin{array}{r}\{\begin{array}{l}{\displaystyle \frac{\mathrm{d}{S}_D(t)}{\mathrm{d}t}={\theta }_D-\frac{S_D(t)}{N_D(t)}{\int }_0^{\mathrm{\infty }}{\beta }_T(a)⟨\mathbf{1},{\mathit{i}}_T(t,a)⟩\mathrm{d}a-{\mu }_D{S}_D(t),}\\ {\displaystyle \frac{\mathrm{\partial }{\mathit{i}}_D(t,a)}{\mathrm{\partial }t}+\frac{\mathrm{\partial }{\mathit{i}}_D(t,a)}{\mathrm{\partial }a}=-({\mu }_D+\mathit{\theta }(a)){\mathit{i}}_D(t,a),}\\ {\mathit{i}}_D(t,0)={\displaystyle \frac{S_D(t)+ϵR_D(t)}{N_D(t)}{\int }_0^{\mathrm{\infty }}{\beta }_T(a)⟨\mathbf{1},{\mathit{i}}_T(t,a)⟩\mathrm{d}a{\mathbf{v}}_1}\\ {\displaystyle +{\int }_0^{\mathrm{\infty }}{\nu }_D(a)⟨{\mathbf{v}}_1,{\mathit{i}}_D(t,a)⟩\mathrm{d}a{\mathbf{v}}_2,}\\ {\displaystyle \frac{\mathrm{d}{R}_D(t)}{\mathrm{d}t}={\int }_0^{\mathrm{\infty }}\gamma (a)⟨{\mathbf{v}}_1,{\mathit{i}}_D(t,a)⟩\mathrm{d}a-ϵ\frac{R_D(t)}{N_D(t)}{\int }_0^{\mathrm{\infty }}{\beta }_T(a)⟨\mathbf{1},{\mathit{i}}_T(t,a)⟩\mathrm{d}a-{\mu }_D{R}_D(t),}\end{array}\end{array}$(7) )–(8(8) $\begin{array}{r}\{\begin{array}{l}{\displaystyle \frac{\mathrm{d}{S}_{TE}(t)}{\mathrm{d}t}=r_E(1-\frac{S_{TE}(t)}{K}){\displaystyle (⟨{\mathit{S}}_T(t),{\mathbf{e}}_3⟩+{\int }_0^{\mathrm{\infty }}⟨{\mathit{i}}_T(t,a),{\mathbf{e}}_3⟩\mathrm{d}a)-({\alpha }_E+{\mu }_E)S_{TE}(t),}}\\ {\displaystyle \frac{\mathrm{d}{\mathit{S}}_T(t)}{\mathrm{d}t}={\alpha }_E{S}_{TE}(t){\mathbf{e}}_1-\mathbf{M}{\mathit{S}}_T(t)-\frac{{\mathit{S}}_T(t)}{N_D(t)}{\int }_0^{\mathrm{\infty }}⟨{\mathit{\beta }}_D(a),{\mathit{i}}_D(t,a)⟩\mathrm{d}a,}\\ {\displaystyle \frac{\mathrm{\partial }{\mathit{i}}_T(t,a)}{\mathrm{\partial }t}+\frac{\mathrm{\partial }{\mathit{i}}_T(t,a)}{\mathrm{\partial }a}=-\mathbf{M}{\mathit{i}}_T(t,a),}\\ {\mathit{i}}_T(t,0)={\displaystyle \frac{{\mathit{S}}_T(t)}{N_D(t)}{\int }_0^{\mathrm{\infty }}⟨{\mathit{\beta }}_D(a),{\mathit{i}}_D(t,a)⟩\mathrm{d}a,}\end{array}\end{array}$(8) ). Moreover, if we set (9) $\begin{array}{r}\mathbf{\Pi }(a,\tau ):=e^{-{\int }_{\tau }^a({\mu }_D+\mathit{\theta }(l))\mathrm{d}l},\mathrm{\forall }a\ge \tau \ge 0.\end{array}$(9) then the solution satisfies, for all $t\ge 0$ and $a\ge 0$, $\begin{array}{rl}& \{\begin{array}{l}{\displaystyle \frac{\mathrm{d}{S}_D(t)}{\mathrm{d}t}={\theta }_D-\frac{S_D(t)}{N_D(t)}{\int }_0^{\mathrm{\infty }}{\beta }_T(a)⟨\mathbf{1},{\mathit{i}}_T(t,a)⟩\mathrm{d}a-{\mu }_D{S}_D(t),}\\ {\displaystyle \frac{\mathrm{d}{R}_D(t)}{\mathrm{d}t}={\int }_0^{\mathrm{\infty }}\gamma (a)i_D(t,a)\mathrm{d}a-ϵ\frac{R_D(t)}{N_D(t)}{\int }_0^{\mathrm{\infty }}{\beta }_T(a)⟨\mathbf{1},{\mathit{i}}_T(t,a)⟩\mathrm{d}a-{\mu }_D{R}_D(t),}\\ {\displaystyle \frac{\mathrm{d}{S}_{TE}(t)}{\mathrm{d}t}=r_E(1-\frac{S_{TE}(t)}{K}){\displaystyle (⟨{\mathit{S}}_T(t),{\mathbf{e}}_3⟩+{\int }_0^{\mathrm{\infty }}⟨{\mathit{i}}_T(t,a),{\mathbf{e}}_3⟩\mathrm{d}a)-({\alpha }_E+{\mu }_E)S_{TE}(t),}}\\ {\displaystyle \frac{\mathrm{d}{\mathit{S}}_T(t)}{\mathrm{d}t}={\alpha }_E{S}_{TE}(t){\mathbf{e}}_1-\mathbf{M}{\mathit{S}}_T(t)-\frac{{\mathit{S}}_T(t)}{N_D(t)}{\int }_0^{\mathrm{\infty }}⟨{\mathit{\beta }}_D(a),{\mathit{i}}_D(t,a)⟩\mathrm{d}a,}\end{array}\\ & {\mathit{i}}_D(t,a)=\{\begin{array}{l}\mathbf{\Pi }(a,a-t){\mathit{i}}_{D0}(a-t),\mathrm{if}a\ge t\\ {\displaystyle \frac{S_D(t-a)+ϵR_D(t-a)}{N_D(t-a)}{\int }_0^{\mathrm{\infty }}{\beta }_T(a)⟨\mathbf{1},{\mathit{i}}_T(t-a,\tau )⟩\mathrm{d}\tau \mathbf{\Pi }(a,0){\mathbf{v}}_1}\\ {\displaystyle +{\int }_0^{\mathrm{\infty }}{\nu }_D(\tau )⟨{\mathbf{v}}_1,{\mathit{i}}_D(t-a,\tau )⟩\mathrm{d}\tau \mathbf{\Pi }(a,0){\mathbf{v}}_2,\mathrm{if}a<t}\end{array}\end{array}$and $\begin{array}{r}{\mathit{i}}_T(t,a)=\{\begin{array}{ll}{e}^{-\mathbf{M}t}{\mathit{i}}_{T0}(a-t)& \mathrm{if}a\ge t\\ e^{-\mathbf{M}a}{\displaystyle \frac{{\mathit{S}}_T(t-a)}{N_D(t-a)}{\int }_0^{\mathrm{\infty }}⟨{\mathit{\beta }}_D(\tau ),{\mathit{i}}_D(t-a,\tau )⟩\mathrm{d}\tau }& \mathrm{if}a<t.\end{array}\end{array}$Moreover, the total population of Dogs, $N_D(t)$, at time t is the unique solution to (10) $\begin{array}{r}\{\begin{array}{l}{\displaystyle \frac{\mathrm{d}{N}_D(t)}{\mathrm{d}t}={\theta }_D-{\mu }_D{N}_D(t)-{\int }_0^{\mathrm{\infty }}{\delta }_D(a)i_D(t,a)\mathrm{d}a-{\int }_0^{\mathrm{\infty }}{\delta }_C(a)c_D(t,a)\mathrm{d}a,t>0}\\ {\displaystyle N_D(0)=S_{TE0}+S_{TL0}+S_{TN0}+S_{TA0}+{\int }_0^{\mathrm{\infty }}[i_{TL0}(a)+i_{TN0}(a)+i_{TA0}(a)]\mathrm{d}a.}\end{array}\end{array}$(10)

Proof.

The existence and positiveness of solutions of system (3(3) $\begin{array}{r}\{\begin{array}{l}{\displaystyle \frac{\mathrm{d}{S}_D(t)}{\mathrm{d}t}={\theta }_D-\frac{S_D(t)}{N_D(t)}{\int }_0^{\mathrm{\infty }}{\lambda }_T(t,a)\mathrm{d}a-{\mu }_D{S}_D(t),}\\ i_D(t,0)={\displaystyle \frac{S_D(t)+ϵR_D(t)}{N_D(t)}{\int }_0^{\mathrm{\infty }}{\lambda }_T(t,a)\mathrm{d}a,}\\ {\displaystyle \frac{\mathrm{\partial }{i}_D(t,a)}{\mathrm{\partial }t}+\frac{\mathrm{\partial }{i}_D(t,a)}{\mathrm{\partial }a}=-({\mu }_D+\gamma (a)+{\delta }_D(a)+{\nu }_D(a))i_D(t,a),}\\ {\displaystyle \frac{\mathrm{\partial }{c}_D(t,a)}{\mathrm{\partial }t}+\frac{\mathrm{\partial }{c}_D(t,a)}{\mathrm{\partial }a}=-({\mu }_D+{\delta }_C(a))c_D(t,a)}\\ c_D(t,0)={\displaystyle {\int }_0^{\mathrm{\infty }}{\nu }_D(a)i_D(t,a)\mathrm{d}a,}\\ {\displaystyle \frac{\mathrm{d}{R}_D(t)}{\mathrm{d}t}={\int }_0^{\mathrm{\infty }}\gamma (a)i_D(t,a)\mathrm{d}a-ϵ\frac{R_D(t)}{N_D(t)}{\int }_0^{\mathrm{\infty }}{\lambda }_T(t,a)\mathrm{d}a-{\mu }_D{R}_D(t).}\end{array}\end{array}$(3) )–(4(4) $\begin{array}{r}\{\begin{array}{l}{\displaystyle \frac{\mathrm{d}{S}_{TE}(t)}{\mathrm{d}t}=r_E(1-\frac{S_{TE}(t)}{K})N_{TA}(t)-({\alpha }_E+{\mu }_E)S_{TE}(t),}\\ {\displaystyle \frac{\mathrm{d}{S}_{TL}(t)}{\mathrm{d}t}={\alpha }_E{S}_{TE}(t)-\frac{S_{TL}(t)}{N_D(t)}{\lambda }_D(t)-({\alpha }_L+{\mu }_L)S_{TL}(t),}\\ {\displaystyle \frac{\mathrm{d}{S}_{TN}(t)}{\mathrm{d}t}={\alpha }_L{S}_{TL}(t)-\frac{S_{TN}(t)}{N_D(t)}{\lambda }_D(t)-({\alpha }_N+{\mu }_N)S_{TN}(t),}\\ {\displaystyle \frac{\mathrm{d}{S}_{TA}(t)}{\mathrm{d}t}={\alpha }_N{S}_{TN}(t)-\frac{S_{TA}(t)}{N_D(t)}{\lambda }_D(t)-{\mu }_A{S}_{TA}(t),}\\ \text{and for }k\in \mathcal{S}\setminus \{E\},\\ i_{Tk}(t,0)={\displaystyle \frac{S_{Tk}(t)}{N_D(t)}{\lambda }_D(t),}\\ {\displaystyle \frac{\mathrm{\partial }{i}_{Tk}(t,a)}{\mathrm{\partial }t}+\frac{\mathrm{\partial }{i}_{Tk}(t,a)}{\mathrm{\partial }a}=-({\mu }_k+{\alpha }_k)i_{Tk}(t,a)}\end{array}\end{array}$(4) ) can be addressed using an integrated semigroup approach and Volterra integral formulation. More precisely, a similar approach as in [37] can be applied for detailed proof of Theorem 3.2. However, here we give a brief sketch of such an approach. For all $t\in \mathbb{R}$, let us set $\begin{array}{r}u(t)=(S_D(t),R_D(t),0_{{\mathbb{R}}^2},{\mathit{i}}_D(t,\cdot ),S_{TE}(t),{\mathit{S}}_T(t),0_{{\mathbb{R}}^3},{\mathit{i}}_T(t,\cdot ))\in \mathcal{Y},\end{array}$where $\mathcal{Y}$ is the space $\begin{array}{r}\mathcal{Y}=\mathbb{R}\times \mathbb{R}\times {\mathbb{R}}^2\times L^1((0,\mathrm{\infty }),{\mathbb{R}}^2)\times {\mathbb{R}}^3\times {\mathbb{R}}^3\times {\mathbb{R}}^3\times L^1((0,\mathrm{\infty }),{\mathbb{R}}^3).\end{array}$Such a Banach space $\mathcal{Y}$ is endow with the usual product norm $\Vert \cdot {\Vert }_{\mathcal{Y}}$.

Let $\mathcal{M}$ : $D(\mathcal{M})\subset \mathcal{Y}⟶\mathcal{Y}$ the linear operator defined by $\begin{array}{r}D(\mathcal{M})=\mathbb{R}\times \mathbb{R}\times \{0_{{\mathbb{R}}^2}\}\times W^{1,1}((0,\mathrm{\infty }),{\mathbb{R}}^2)\times {\mathbb{R}}^3\times {\mathbb{R}}^3\times \{0_{{\mathbb{R}}^3}\}\times W^{1,1}((0,\mathrm{\infty }),{\mathbb{R}}^3),\end{array}$and $\begin{array}{r}\mathcal{M}\left(\begin{array}{c}{S}_D\\ R_D\\ 0_{{\mathbb{R}}^2}\\ {\mathit{i}}_D(\cdot )\\ S_{TE}\\ {\mathit{S}}_T\\ 0_{{\mathbb{R}}^3}\\ {\mathit{i}}_T(\cdot )\end{array}\right)=\left(\begin{array}{c}-{\mu }_D{S}_D\\ -{\mu }_D{R}_D\\ -{\mathit{i}}_D(0)\\ -{\mathit{i}}_D^{\prime }-({\mu }_D+\mathit{\theta }){\mathit{i}}_D\\ -({\alpha }_E+{\mu }_E)S_{TE}\\ -\mathbf{M}{\mathit{S}}_T\\ -{\mathit{i}}_T(0)\\ -{\mathit{i}}_T^{\prime }-\mathbf{M}{\mathit{i}}_T\end{array}\right),\end{array}$as well as the map $\mathcal{F}:\mathcal{Y}\to \mathcal{Y}$ such that $\begin{array}{r}\mathcal{F}\left(\begin{array}{c}{S}_D\\ R_D\\ 0_{{\mathbb{R}}^2}\\ {\mathit{i}}_D(\cdot )\\ S_{TE}\\ {\mathit{S}}_T\\ 0_{{\mathbb{R}}^3}\\ {\mathit{i}}_T(\cdot )\end{array}\right)=\left(\begin{array}{c}{\displaystyle {\theta }_D-\frac{S_D}{N_D}{\int }_0^{\mathrm{\infty }}{\beta }_T(a)⟨\mathbf{1},{\mathit{i}}_T(a)⟩\mathrm{d}a}\\ {\displaystyle {\int }_0^{\mathrm{\infty }}\gamma (a)⟨{\mathbf{v}}_1,{\mathit{i}}_D(a)⟩\mathrm{d}a-ϵ\frac{R_D}{N_D}{\int }_0^{\mathrm{\infty }}{\beta }_T(a)⟨\mathbf{1},{\mathit{i}}_T(a)⟩\mathrm{d}a}\\ {\displaystyle \frac{S_D+ϵR_D}{N_D}{\int }_0^{\mathrm{\infty }}{\beta }_T(a)⟨\mathbf{1},{\mathit{i}}_T(a)⟩\mathrm{d}a{\mathbf{v}}_1+{\int }_0^{\mathrm{\infty }}{\nu }_D(a)⟨{\mathbf{v}}_1,{\mathit{i}}_D(a)⟩\mathrm{d}a{\mathbf{v}}_2}\\ 0\\ {\displaystyle r_E(1-\frac{S_{TE}}{K})(⟨{\mathit{S}}_T,{\mathbf{e}}_3⟩+{\int }_0^{\mathrm{\infty }}⟨{\mathit{i}}_T(a),{\mathbf{e}}_3⟩\mathrm{d}a)}\\ {\displaystyle {\alpha }_E{S}_{TE}{\mathbf{e}}_1-\frac{{\mathit{S}}_T}{N_D}{\int }_0^{\mathrm{\infty }}⟨{\mathit{\beta }}_D(a),{\mathit{i}}_D(a)⟩\mathrm{d}a}\\ {\displaystyle \frac{{\mathit{S}}_T}{N_D}{\int }_0^{\mathrm{\infty }}⟨{\mathit{\beta }}_D(a),{\mathit{i}}_D(a)⟩\mathrm{d}a}\\ 0\end{array}\right).\end{array}$Observe that the nonlinear map $\mathcal{F}$ is not well defined on $\overline{D(\mathcal{M})}$ due to the term $N_D$ and is then not locally Lipschitz continuous. Furthermore, for any $\zeta >0$, let us introduce the space $\begin{array}{r}{\mathcal{Y}}_{\zeta }=\{u\in \overline{D(\mathcal{M})}:\kappa (u)\ge \zeta \}\subset \overline{D(\mathcal{M})},\end{array}$where $\kappa :\mathcal{Y}\to \mathcal{Y}$ is the operator defined by $\begin{array}{r}\kappa (u)=S_D+{\int }_0^{\mathrm{\infty }}{i}_D(a)\mathrm{d}a+{\int }_0^{\mathrm{\infty }}{c}_D(a)\mathrm{d}a+R_D.\end{array}$We can now define the nonlinear operator $F_{\zeta }:{\mathcal{Y}}_{\zeta }\to \mathcal{Y}$ by $F_{\zeta }\equiv F$. Therefore, system (7(7) $\begin{array}{r}\{\begin{array}{l}{\displaystyle \frac{\mathrm{d}{S}_D(t)}{\mathrm{d}t}={\theta }_D-\frac{S_D(t)}{N_D(t)}{\int }_0^{\mathrm{\infty }}{\beta }_T(a)⟨\mathbf{1},{\mathit{i}}_T(t,a)⟩\mathrm{d}a-{\mu }_D{S}_D(t),}\\ {\displaystyle \frac{\mathrm{\partial }{\mathit{i}}_D(t,a)}{\mathrm{\partial }t}+\frac{\mathrm{\partial }{\mathit{i}}_D(t,a)}{\mathrm{\partial }a}=-({\mu }_D+\mathit{\theta }(a)){\mathit{i}}_D(t,a),}\\ {\mathit{i}}_D(t,0)={\displaystyle \frac{S_D(t)+ϵR_D(t)}{N_D(t)}{\int }_0^{\mathrm{\infty }}{\beta }_T(a)⟨\mathbf{1},{\mathit{i}}_T(t,a)⟩\mathrm{d}a{\mathbf{v}}_1}\\ {\displaystyle +{\int }_0^{\mathrm{\infty }}{\nu }_D(a)⟨{\mathbf{v}}_1,{\mathit{i}}_D(t,a)⟩\mathrm{d}a{\mathbf{v}}_2,}\\ {\displaystyle \frac{\mathrm{d}{R}_D(t)}{\mathrm{d}t}={\int }_0^{\mathrm{\infty }}\gamma (a)⟨{\mathbf{v}}_1,{\mathit{i}}_D(t,a)⟩\mathrm{d}a-ϵ\frac{R_D(t)}{N_D(t)}{\int }_0^{\mathrm{\infty }}{\beta }_T(a)⟨\mathbf{1},{\mathit{i}}_T(t,a)⟩\mathrm{d}a-{\mu }_D{R}_D(t),}\end{array}\end{array}$(7) )–(8(8) $\begin{array}{r}\{\begin{array}{l}{\displaystyle \frac{\mathrm{d}{S}_{TE}(t)}{\mathrm{d}t}=r_E(1-\frac{S_{TE}(t)}{K}){\displaystyle (⟨{\mathit{S}}_T(t),{\mathbf{e}}_3⟩+{\int }_0^{\mathrm{\infty }}⟨{\mathit{i}}_T(t,a),{\mathbf{e}}_3⟩\mathrm{d}a)-({\alpha }_E+{\mu }_E)S_{TE}(t),}}\\ {\displaystyle \frac{\mathrm{d}{\mathit{S}}_T(t)}{\mathrm{d}t}={\alpha }_E{S}_{TE}(t){\mathbf{e}}_1-\mathbf{M}{\mathit{S}}_T(t)-\frac{{\mathit{S}}_T(t)}{N_D(t)}{\int }_0^{\mathrm{\infty }}⟨{\mathit{\beta }}_D(a),{\mathit{i}}_D(t,a)⟩\mathrm{d}a,}\\ {\displaystyle \frac{\mathrm{\partial }{\mathit{i}}_T(t,a)}{\mathrm{\partial }t}+\frac{\mathrm{\partial }{\mathit{i}}_T(t,a)}{\mathrm{\partial }a}=-\mathbf{M}{\mathit{i}}_T(t,a),}\\ {\mathit{i}}_T(t,0)={\displaystyle \frac{{\mathit{S}}_T(t)}{N_D(t)}{\int }_0^{\mathrm{\infty }}⟨{\mathit{\beta }}_D(a),{\mathit{i}}_D(t,a)⟩\mathrm{d}a,}\end{array}\end{array}$(8) ) can be rewritten as the following non-densely defined abstract Cauchy problem: (11) $\begin{array}{r}\frac{\mathrm{d}}{\mathrm{d}t}u(t)=\mathcal{M}u(t)+{\mathcal{F}}_{\zeta }(u(t)),t>0,u(0)\in {\mathcal{Y}}_{\zeta }\cap {\mathcal{Y}}_{+}.\end{array}$(11) We then address the existence and uniqueness of bounded solutions to (11(11) $\begin{array}{r}\frac{\mathrm{d}}{\mathrm{d}t}u(t)=\mathcal{M}u(t)+{\mathcal{F}}_{\zeta }(u(t)),t>0,u(0)\in {\mathcal{Y}}_{\zeta }\cap {\mathcal{Y}}_{+}.\end{array}$(11) ) in a manner similar to the approach in [37].

3.2. Disease invasion process

System (7(7) $\begin{array}{r}\{\begin{array}{l}{\displaystyle \frac{\mathrm{d}{S}_D(t)}{\mathrm{d}t}={\theta }_D-\frac{S_D(t)}{N_D(t)}{\int }_0^{\mathrm{\infty }}{\beta }_T(a)⟨\mathbf{1},{\mathit{i}}_T(t,a)⟩\mathrm{d}a-{\mu }_D{S}_D(t),}\\ {\displaystyle \frac{\mathrm{\partial }{\mathit{i}}_D(t,a)}{\mathrm{\partial }t}+\frac{\mathrm{\partial }{\mathit{i}}_D(t,a)}{\mathrm{\partial }a}=-({\mu }_D+\mathit{\theta }(a)){\mathit{i}}_D(t,a),}\\ {\mathit{i}}_D(t,0)={\displaystyle \frac{S_D(t)+ϵR_D(t)}{N_D(t)}{\int }_0^{\mathrm{\infty }}{\beta }_T(a)⟨\mathbf{1},{\mathit{i}}_T(t,a)⟩\mathrm{d}a{\mathbf{v}}_1}\\ {\displaystyle +{\int }_0^{\mathrm{\infty }}{\nu }_D(a)⟨{\mathbf{v}}_1,{\mathit{i}}_D(t,a)⟩\mathrm{d}a{\mathbf{v}}_2,}\\ {\displaystyle \frac{\mathrm{d}{R}_D(t)}{\mathrm{d}t}={\int }_0^{\mathrm{\infty }}\gamma (a)⟨{\mathbf{v}}_1,{\mathit{i}}_D(t,a)⟩\mathrm{d}a-ϵ\frac{R_D(t)}{N_D(t)}{\int }_0^{\mathrm{\infty }}{\beta }_T(a)⟨\mathbf{1},{\mathit{i}}_T(t,a)⟩\mathrm{d}a-{\mu }_D{R}_D(t),}\end{array}\end{array}$(7) )–(8(8) $\begin{array}{r}\{\begin{array}{l}{\displaystyle \frac{\mathrm{d}{S}_{TE}(t)}{\mathrm{d}t}=r_E(1-\frac{S_{TE}(t)}{K}){\displaystyle (⟨{\mathit{S}}_T(t),{\mathbf{e}}_3⟩+{\int }_0^{\mathrm{\infty }}⟨{\mathit{i}}_T(t,a),{\mathbf{e}}_3⟩\mathrm{d}a)-({\alpha }_E+{\mu }_E)S_{TE}(t),}}\\ {\displaystyle \frac{\mathrm{d}{\mathit{S}}_T(t)}{\mathrm{d}t}={\alpha }_E{S}_{TE}(t){\mathbf{e}}_1-\mathbf{M}{\mathit{S}}_T(t)-\frac{{\mathit{S}}_T(t)}{N_D(t)}{\int }_0^{\mathrm{\infty }}⟨{\mathit{\beta }}_D(a),{\mathit{i}}_D(t,a)⟩\mathrm{d}a,}\\ {\displaystyle \frac{\mathrm{\partial }{\mathit{i}}_T(t,a)}{\mathrm{\partial }t}+\frac{\mathrm{\partial }{\mathit{i}}_T(t,a)}{\mathrm{\partial }a}=-\mathbf{M}{\mathit{i}}_T(t,a),}\\ {\mathit{i}}_T(t,0)={\displaystyle \frac{{\mathit{S}}_T(t)}{N_D(t)}{\int }_0^{\mathrm{\infty }}⟨{\mathit{\beta }}_D(a),{\mathit{i}}_D(t,a)⟩\mathrm{d}a,}\end{array}\end{array}$(8) ) exhibits two disease-free equilibria. More precisely, in an infection-free environment, we have a disease-free equilibrium either in the presence of dogs and ticks or only in the presence of dogs. But, only the former is interesting within this context because the disease is transmitted from dogs to ticks and vice versa. Such an equilibrium is named after the disease-free equilibrium of system (3(3) $\begin{array}{r}\{\begin{array}{l}{\displaystyle \frac{\mathrm{d}{S}_D(t)}{\mathrm{d}t}={\theta }_D-\frac{S_D(t)}{N_D(t)}{\int }_0^{\mathrm{\infty }}{\lambda }_T(t,a)\mathrm{d}a-{\mu }_D{S}_D(t),}\\ i_D(t,0)={\displaystyle \frac{S_D(t)+ϵR_D(t)}{N_D(t)}{\int }_0^{\mathrm{\infty }}{\lambda }_T(t,a)\mathrm{d}a,}\\ {\displaystyle \frac{\mathrm{\partial }{i}_D(t,a)}{\mathrm{\partial }t}+\frac{\mathrm{\partial }{i}_D(t,a)}{\mathrm{\partial }a}=-({\mu }_D+\gamma (a)+{\delta }_D(a)+{\nu }_D(a))i_D(t,a),}\\ {\displaystyle \frac{\mathrm{\partial }{c}_D(t,a)}{\mathrm{\partial }t}+\frac{\mathrm{\partial }{c}_D(t,a)}{\mathrm{\partial }a}=-({\mu }_D+{\delta }_C(a))c_D(t,a)}\\ c_D(t,0)={\displaystyle {\int }_0^{\mathrm{\infty }}{\nu }_D(a)i_D(t,a)\mathrm{d}a,}\\ {\displaystyle \frac{\mathrm{d}{R}_D(t)}{\mathrm{d}t}={\int }_0^{\mathrm{\infty }}\gamma (a)i_D(t,a)\mathrm{d}a-ϵ\frac{R_D(t)}{N_D(t)}{\int }_0^{\mathrm{\infty }}{\lambda }_T(t,a)\mathrm{d}a-{\mu }_D{R}_D(t).}\end{array}\end{array}$(3) )–(4(4) $\begin{array}{r}\{\begin{array}{l}{\displaystyle \frac{\mathrm{d}{S}_{TE}(t)}{\mathrm{d}t}=r_E(1-\frac{S_{TE}(t)}{K})N_{TA}(t)-({\alpha }_E+{\mu }_E)S_{TE}(t),}\\ {\displaystyle \frac{\mathrm{d}{S}_{TL}(t)}{\mathrm{d}t}={\alpha }_E{S}_{TE}(t)-\frac{S_{TL}(t)}{N_D(t)}{\lambda }_D(t)-({\alpha }_L+{\mu }_L)S_{TL}(t),}\\ {\displaystyle \frac{\mathrm{d}{S}_{TN}(t)}{\mathrm{d}t}={\alpha }_L{S}_{TL}(t)-\frac{S_{TN}(t)}{N_D(t)}{\lambda }_D(t)-({\alpha }_N+{\mu }_N)S_{TN}(t),}\\ {\displaystyle \frac{\mathrm{d}{S}_{TA}(t)}{\mathrm{d}t}={\alpha }_N{S}_{TN}(t)-\frac{S_{TA}(t)}{N_D(t)}{\lambda }_D(t)-{\mu }_A{S}_{TA}(t),}\\ \text{and for }k\in \mathcal{S}\setminus \{E\},\\ i_{Tk}(t,0)={\displaystyle \frac{S_{Tk}(t)}{N_D(t)}{\lambda }_D(t),}\\ {\displaystyle \frac{\mathrm{\partial }{i}_{Tk}(t,a)}{\mathrm{\partial }t}+\frac{\mathrm{\partial }{i}_{Tk}(t,a)}{\mathrm{\partial }a}=-({\mu }_k+{\alpha }_k)i_{Tk}(t,a)}\end{array}\end{array}$(4) ). Let us introduce the following threshold parameter (12) $\begin{array}{r}{\mathcal{R}}_{ST}^0=\frac{r_E}{{\mu }_A}\prod _{k\in \mathcal{S}}\frac{{\alpha }_k}{{\alpha }_k+{\mu }_k}=r_E\frac{{\alpha }_E}{{\alpha }_E+{\mu }_E}\frac{{\alpha }_L}{{\alpha }_L+{\mu }_L}\frac{{\alpha }_N}{{\alpha }_N+{\mu }_N}\frac{1}{{\mu }_A}.\end{array}$(12) The parameter ${\mathcal{R}}_{ST}^0$ accounts for the ticks' reproduction number and is such that, $\frac{{\alpha }_E}{{\alpha }_E+{\mu }_E}$ represents the fraction of eggs that progress to the larval stage, $\frac{{\alpha }_L}{{\alpha }_L+{\mu }_L}$ is the fraction of larval that progress to the nymphal stage, $\frac{{\alpha }_N}{{\alpha }_N+{\mu }_N}$ is the fraction of nymph that progress to the adult stage, $\frac{1}{{\mu }_A}$ is the life expectancy of adult ticks and $r_E$ is the ticks egg laying rate. Furthermore, when ${\mathcal{R}}_{ST}^0>1$, the disease-free equilibrium of system (3(3) $\begin{array}{r}\{\begin{array}{l}{\displaystyle \frac{\mathrm{d}{S}_D(t)}{\mathrm{d}t}={\theta }_D-\frac{S_D(t)}{N_D(t)}{\int }_0^{\mathrm{\infty }}{\lambda }_T(t,a)\mathrm{d}a-{\mu }_D{S}_D(t),}\\ i_D(t,0)={\displaystyle \frac{S_D(t)+ϵR_D(t)}{N_D(t)}{\int }_0^{\mathrm{\infty }}{\lambda }_T(t,a)\mathrm{d}a,}\\ {\displaystyle \frac{\mathrm{\partial }{i}_D(t,a)}{\mathrm{\partial }t}+\frac{\mathrm{\partial }{i}_D(t,a)}{\mathrm{\partial }a}=-({\mu }_D+\gamma (a)+{\delta }_D(a)+{\nu }_D(a))i_D(t,a),}\\ {\displaystyle \frac{\mathrm{\partial }{c}_D(t,a)}{\mathrm{\partial }t}+\frac{\mathrm{\partial }{c}_D(t,a)}{\mathrm{\partial }a}=-({\mu }_D+{\delta }_C(a))c_D(t,a)}\\ c_D(t,0)={\displaystyle {\int }_0^{\mathrm{\infty }}{\nu }_D(a)i_D(t,a)\mathrm{d}a,}\\ {\displaystyle \frac{\mathrm{d}{R}_D(t)}{\mathrm{d}t}={\int }_0^{\mathrm{\infty }}\gamma (a)i_D(t,a)\mathrm{d}a-ϵ\frac{R_D(t)}{N_D(t)}{\int }_0^{\mathrm{\infty }}{\lambda }_T(t,a)\mathrm{d}a-{\mu }_D{R}_D(t).}\end{array}\end{array}$(3) )–(4(4) $\begin{array}{r}\{\begin{array}{l}{\displaystyle \frac{\mathrm{d}{S}_{TE}(t)}{\mathrm{d}t}=r_E(1-\frac{S_{TE}(t)}{K})N_{TA}(t)-({\alpha }_E+{\mu }_E)S_{TE}(t),}\\ {\displaystyle \frac{\mathrm{d}{S}_{TL}(t)}{\mathrm{d}t}={\alpha }_E{S}_{TE}(t)-\frac{S_{TL}(t)}{N_D(t)}{\lambda }_D(t)-({\alpha }_L+{\mu }_L)S_{TL}(t),}\\ {\displaystyle \frac{\mathrm{d}{S}_{TN}(t)}{\mathrm{d}t}={\alpha }_L{S}_{TL}(t)-\frac{S_{TN}(t)}{N_D(t)}{\lambda }_D(t)-({\alpha }_N+{\mu }_N)S_{TN}(t),}\\ {\displaystyle \frac{\mathrm{d}{S}_{TA}(t)}{\mathrm{d}t}={\alpha }_N{S}_{TN}(t)-\frac{S_{TA}(t)}{N_D(t)}{\lambda }_D(t)-{\mu }_A{S}_{TA}(t),}\\ \text{and for }k\in \mathcal{S}\setminus \{E\},\\ i_{Tk}(t,0)={\displaystyle \frac{S_{Tk}(t)}{N_D(t)}{\lambda }_D(t),}\\ {\displaystyle \frac{\mathrm{\partial }{i}_{Tk}(t,a)}{\mathrm{\partial }t}+\frac{\mathrm{\partial }{i}_{Tk}(t,a)}{\mathrm{\partial }a}=-({\mu }_k+{\alpha }_k)i_{Tk}(t,a)}\end{array}\end{array}$(4) ) is given by (13) $\begin{array}{r}{\mathcal{E}}_0:=({\overline{S}}_D,0_{L^1((0,+\mathrm{\infty }),{\mathbb{R}}^2)},0,{\overline{S}}_{TE},{\overline{\mathit{S}}}_T,0_{L^1((0,+\mathrm{\infty }),{\mathbb{R}}^3)}),\end{array}$(13) where (14) $\begin{array}{r}\{\begin{array}{l}{\displaystyle {\overline{S}}_D=\frac{{\theta }_D}{{\mu }_D},}\\ {\displaystyle {\overline{S}}_{TE}=K(1-\frac{1}{{\mathcal{R}}_{ST}^0}),}\\ {\displaystyle {\overline{\mathit{S}}}_T={\alpha }_E{\overline{S}}_{TE}\sum _{k=1}^3{\mathbf{e}}_k\prod _{j=1}^k\frac{{\alpha }_{j-1}}{{\mu }_j+{\alpha }_j},\mathrm{with}{\alpha }_0=1,{\alpha }_3=0.}\end{array}\end{array}$(14) In the last equality, the correspondence $L\leftarrow 1$, $N\leftarrow 2$, and $A\leftarrow 3$ is used to facilitate the notations. See Section 6 for a detailed computation of ${\mathcal{E}}_0$.

Next, we introduce the following threshold (15) $\begin{array}{r}{\mathcal{T}}_0={\mathcal{T}}_0^{D\to T}\times {\mathcal{T}}_0^{T\to D},\end{array}$(15) where ${\mathcal{T}}_0^{D\to T}$ quantifies the transmission capability from dogs to ticks and ${\mathcal{T}}_0^{T\to D}$ the transmission capability from ticks to dogs. More precisely, we have $\begin{array}{r}{\mathcal{T}}_0^{D\to T}=\underset{\begin{array}{c}\text{Transmission from infected dogs to ticks}\end{array}}{\underbrace{{\int }_0^{\mathrm{\infty }}{\beta }_D^1(a){\pi }_D(a)\mathrm{d}a}}+\underset{\begin{array}{c}\text{Transmission from chronically infected dogs to ticks}\end{array}}{\underbrace{{\int }_0^{\mathrm{\infty }}{\beta }_D^2(a){\pi }_C(a)\mathrm{d}a{\int }_0^{\mathrm{\infty }}{\nu }_D(a){\pi }_D(a)\mathrm{d}a}},\end{array}$and (16) $\begin{array}{r}{\mathcal{T}}_0^{T\to D}=\sum _{k\in \mathcal{S}\setminus \{E\}}\underset{\begin{array}{c}\text{Proportion ticks of stage }\$k\$\end{array}}{\underbrace{\frac{{\overline{S}}_k}{\sum _{j\in \mathcal{S}\setminus \{E\}}{\overline{S}}_j}}}\times {\mathcal{T}}_{0k}^{T\to D},\end{array}$(16) where ${\mathcal{T}}_{0k}^{T\to D}$ denotes the vectorial capacity of ticks population of stage k and is explicitly given by (17) $\begin{array}{r}{\mathcal{T}}_{0k}^{T\to D}=\underset{\begin{array}{c}\text{Ticks/Dogs ratio}\end{array}}{\underbrace{\frac{\sum _{j\in \mathcal{S}\setminus \{E\}}{\overline{S}}_j}{{\overline{S}}_D}}}\times \underset{\begin{array}{c}\text{Transmission from ticks at stage }k\text{ to Dogs}\end{array}}{\underbrace{{\int }_0^{\mathrm{\infty }}{\beta }_T(a)⟨\mathbf{1},e^{-\mathbf{M}a}{\mathbf{e}}_k⟩\mathrm{d}a,}}\end{array}$(17) with (18) $\begin{array}{r}\{\begin{array}{l}{\pi }_D(a)=e^{-{\int }_0^a({\mu }_D+{\delta }_D(l)+{\nu }_D(l))\mathrm{d}l}=\mathbf{\Pi }(a,0){\mathbf{v}}_1,\\ {\pi }_C(a)=e^{-{\int }_0^a({\mu }_D+{\delta }_C(l))\mathrm{d}l}=\mathbf{\Pi }(a,0){\mathbf{v}}_2.\end{array}\end{array}$(18) The map $a\to {\pi }_D(a)$ (resp. $a\to {\pi }_C(a)$) describes the probability to still be infected (resp. chronically infected) a-time post-infection. Based on the above notations, we now state the invasion dynamics in terms of the threshold ${\mathcal{T}}_0$ defined in (15(15) $\begin{array}{r}{\mathcal{T}}_0={\mathcal{T}}_0^{D\to T}\times {\mathcal{T}}_0^{T\to D},\end{array}$(15) ).

Theorem 3.3

Let Assumption 3.1 be satisfied. Assume in addition that the ticks' reproduction number satisfies ${\mathcal{R}}_{ST}^0>1$. Then the following properties hold true:

(i)	If ${\mathcal{T}}_0<1$ then the disease-free equilibrium ${\mathcal{E}}_0$ is locally asymptotically stable.
(ii)	If ${\mathcal{T}}_0>1$ then the disease-free equilibrium ${\mathcal{E}}_0$ is unstable.

The proof of Theorem 3.3 is given in Section 7.

3.3. Existence of an endemic equilibrium and bifurcation

In this section, we state our main result concerning necessary and sufficient conditions for the existence of an endemic equilibrium to system (7(7) $\begin{array}{r}\{\begin{array}{l}{\displaystyle \frac{\mathrm{d}{S}_D(t)}{\mathrm{d}t}={\theta }_D-\frac{S_D(t)}{N_D(t)}{\int }_0^{\mathrm{\infty }}{\beta }_T(a)⟨\mathbf{1},{\mathit{i}}_T(t,a)⟩\mathrm{d}a-{\mu }_D{S}_D(t),}\\ {\displaystyle \frac{\mathrm{\partial }{\mathit{i}}_D(t,a)}{\mathrm{\partial }t}+\frac{\mathrm{\partial }{\mathit{i}}_D(t,a)}{\mathrm{\partial }a}=-({\mu }_D+\mathit{\theta }(a)){\mathit{i}}_D(t,a),}\\ {\mathit{i}}_D(t,0)={\displaystyle \frac{S_D(t)+ϵR_D(t)}{N_D(t)}{\int }_0^{\mathrm{\infty }}{\beta }_T(a)⟨\mathbf{1},{\mathit{i}}_T(t,a)⟩\mathrm{d}a{\mathbf{v}}_1}\\ {\displaystyle +{\int }_0^{\mathrm{\infty }}{\nu }_D(a)⟨{\mathbf{v}}_1,{\mathit{i}}_D(t,a)⟩\mathrm{d}a{\mathbf{v}}_2,}\\ {\displaystyle \frac{\mathrm{d}{R}_D(t)}{\mathrm{d}t}={\int }_0^{\mathrm{\infty }}\gamma (a)⟨{\mathbf{v}}_1,{\mathit{i}}_D(t,a)⟩\mathrm{d}a-ϵ\frac{R_D(t)}{N_D(t)}{\int }_0^{\mathrm{\infty }}{\beta }_T(a)⟨\mathbf{1},{\mathit{i}}_T(t,a)⟩\mathrm{d}a-{\mu }_D{R}_D(t),}\end{array}\end{array}$(7) )–(8(8) $\begin{array}{r}\{\begin{array}{l}{\displaystyle \frac{\mathrm{d}{S}_{TE}(t)}{\mathrm{d}t}=r_E(1-\frac{S_{TE}(t)}{K}){\displaystyle (⟨{\mathit{S}}_T(t),{\mathbf{e}}_3⟩+{\int }_0^{\mathrm{\infty }}⟨{\mathit{i}}_T(t,a),{\mathbf{e}}_3⟩\mathrm{d}a)-({\alpha }_E+{\mu }_E)S_{TE}(t),}}\\ {\displaystyle \frac{\mathrm{d}{\mathit{S}}_T(t)}{\mathrm{d}t}={\alpha }_E{S}_{TE}(t){\mathbf{e}}_1-\mathbf{M}{\mathit{S}}_T(t)-\frac{{\mathit{S}}_T(t)}{N_D(t)}{\int }_0^{\mathrm{\infty }}⟨{\mathit{\beta }}_D(a),{\mathit{i}}_D(t,a)⟩\mathrm{d}a,}\\ {\displaystyle \frac{\mathrm{\partial }{\mathit{i}}_T(t,a)}{\mathrm{\partial }t}+\frac{\mathrm{\partial }{\mathit{i}}_T(t,a)}{\mathrm{\partial }a}=-\mathbf{M}{\mathit{i}}_T(t,a),}\\ {\mathit{i}}_T(t,0)={\displaystyle \frac{{\mathit{S}}_T(t)}{N_D(t)}{\int }_0^{\mathrm{\infty }}⟨{\mathit{\beta }}_D(a),{\mathit{i}}_D(t,a)⟩\mathrm{d}a,}\end{array}\end{array}$(8) ) and forward (resp. backward) bifurcation at ${\mathcal{T}}_0=1$. Denote by ${\mathcal{E}}^{\ast }$ an endemic equilibrium to system (7(7) $\begin{array}{r}\{\begin{array}{l}{\displaystyle \frac{\mathrm{d}{S}_D(t)}{\mathrm{d}t}={\theta }_D-\frac{S_D(t)}{N_D(t)}{\int }_0^{\mathrm{\infty }}{\beta }_T(a)⟨\mathbf{1},{\mathit{i}}_T(t,a)⟩\mathrm{d}a-{\mu }_D{S}_D(t),}\\ {\displaystyle \frac{\mathrm{\partial }{\mathit{i}}_D(t,a)}{\mathrm{\partial }t}+\frac{\mathrm{\partial }{\mathit{i}}_D(t,a)}{\mathrm{\partial }a}=-({\mu }_D+\mathit{\theta }(a)){\mathit{i}}_D(t,a),}\\ {\mathit{i}}_D(t,0)={\displaystyle \frac{S_D(t)+ϵR_D(t)}{N_D(t)}{\int }_0^{\mathrm{\infty }}{\beta }_T(a)⟨\mathbf{1},{\mathit{i}}_T(t,a)⟩\mathrm{d}a{\mathbf{v}}_1}\\ {\displaystyle +{\int }_0^{\mathrm{\infty }}{\nu }_D(a)⟨{\mathbf{v}}_1,{\mathit{i}}_D(t,a)⟩\mathrm{d}a{\mathbf{v}}_2,}\\ {\displaystyle \frac{\mathrm{d}{R}_D(t)}{\mathrm{d}t}={\int }_0^{\mathrm{\infty }}\gamma (a)⟨{\mathbf{v}}_1,{\mathit{i}}_D(t,a)⟩\mathrm{d}a-ϵ\frac{R_D(t)}{N_D(t)}{\int }_0^{\mathrm{\infty }}{\beta }_T(a)⟨\mathbf{1},{\mathit{i}}_T(t,a)⟩\mathrm{d}a-{\mu }_D{R}_D(t),}\end{array}\end{array}$(7) )–(8(8) $\begin{array}{r}\{\begin{array}{l}{\displaystyle \frac{\mathrm{d}{S}_{TE}(t)}{\mathrm{d}t}=r_E(1-\frac{S_{TE}(t)}{K}){\displaystyle (⟨{\mathit{S}}_T(t),{\mathbf{e}}_3⟩+{\int }_0^{\mathrm{\infty }}⟨{\mathit{i}}_T(t,a),{\mathbf{e}}_3⟩\mathrm{d}a)-({\alpha }_E+{\mu }_E)S_{TE}(t),}}\\ {\displaystyle \frac{\mathrm{d}{\mathit{S}}_T(t)}{\mathrm{d}t}={\alpha }_E{S}_{TE}(t){\mathbf{e}}_1-\mathbf{M}{\mathit{S}}_T(t)-\frac{{\mathit{S}}_T(t)}{N_D(t)}{\int }_0^{\mathrm{\infty }}⟨{\mathit{\beta }}_D(a),{\mathit{i}}_D(t,a)⟩\mathrm{d}a,}\\ {\displaystyle \frac{\mathrm{\partial }{\mathit{i}}_T(t,a)}{\mathrm{\partial }t}+\frac{\mathrm{\partial }{\mathit{i}}_T(t,a)}{\mathrm{\partial }a}=-\mathbf{M}{\mathit{i}}_T(t,a),}\\ {\mathit{i}}_T(t,0)={\displaystyle \frac{{\mathit{S}}_T(t)}{N_D(t)}{\int }_0^{\mathrm{\infty }}⟨{\mathit{\beta }}_D(a),{\mathit{i}}_D(t,a)⟩\mathrm{d}a,}\end{array}\end{array}$(8) ) that is a time-independent solution to system (7(7) $\begin{array}{r}\{\begin{array}{l}{\displaystyle \frac{\mathrm{d}{S}_D(t)}{\mathrm{d}t}={\theta }_D-\frac{S_D(t)}{N_D(t)}{\int }_0^{\mathrm{\infty }}{\beta }_T(a)⟨\mathbf{1},{\mathit{i}}_T(t,a)⟩\mathrm{d}a-{\mu }_D{S}_D(t),}\\ {\displaystyle \frac{\mathrm{\partial }{\mathit{i}}_D(t,a)}{\mathrm{\partial }t}+\frac{\mathrm{\partial }{\mathit{i}}_D(t,a)}{\mathrm{\partial }a}=-({\mu }_D+\mathit{\theta }(a)){\mathit{i}}_D(t,a),}\\ {\mathit{i}}_D(t,0)={\displaystyle \frac{S_D(t)+ϵR_D(t)}{N_D(t)}{\int }_0^{\mathrm{\infty }}{\beta }_T(a)⟨\mathbf{1},{\mathit{i}}_T(t,a)⟩\mathrm{d}a{\mathbf{v}}_1}\\ {\displaystyle +{\int }_0^{\mathrm{\infty }}{\nu }_D(a)⟨{\mathbf{v}}_1,{\mathit{i}}_D(t,a)⟩\mathrm{d}a{\mathbf{v}}_2,}\\ {\displaystyle \frac{\mathrm{d}{R}_D(t)}{\mathrm{d}t}={\int }_0^{\mathrm{\infty }}\gamma (a)⟨{\mathbf{v}}_1,{\mathit{i}}_D(t,a)⟩\mathrm{d}a-ϵ\frac{R_D(t)}{N_D(t)}{\int }_0^{\mathrm{\infty }}{\beta }_T(a)⟨\mathbf{1},{\mathit{i}}_T(t,a)⟩\mathrm{d}a-{\mu }_D{R}_D(t),}\end{array}\end{array}$(7) )–(8(8) $\begin{array}{r}\{\begin{array}{l}{\displaystyle \frac{\mathrm{d}{S}_{TE}(t)}{\mathrm{d}t}=r_E(1-\frac{S_{TE}(t)}{K}){\displaystyle (⟨{\mathit{S}}_T(t),{\mathbf{e}}_3⟩+{\int }_0^{\mathrm{\infty }}⟨{\mathit{i}}_T(t,a),{\mathbf{e}}_3⟩\mathrm{d}a)-({\alpha }_E+{\mu }_E)S_{TE}(t),}}\\ {\displaystyle \frac{\mathrm{d}{\mathit{S}}_T(t)}{\mathrm{d}t}={\alpha }_E{S}_{TE}(t){\mathbf{e}}_1-\mathbf{M}{\mathit{S}}_T(t)-\frac{{\mathit{S}}_T(t)}{N_D(t)}{\int }_0^{\mathrm{\infty }}⟨{\mathit{\beta }}_D(a),{\mathit{i}}_D(t,a)⟩\mathrm{d}a,}\\ {\displaystyle \frac{\mathrm{\partial }{\mathit{i}}_T(t,a)}{\mathrm{\partial }t}+\frac{\mathrm{\partial }{\mathit{i}}_T(t,a)}{\mathrm{\partial }a}=-\mathbf{M}{\mathit{i}}_T(t,a),}\\ {\mathit{i}}_T(t,0)={\displaystyle \frac{{\mathit{S}}_T(t)}{N_D(t)}{\int }_0^{\mathrm{\infty }}⟨{\mathit{\beta }}_D(a),{\mathit{i}}_D(t,a)⟩\mathrm{d}a,}\end{array}\end{array}$(8) ). (19) $\begin{array}{r}{\mathcal{E}}^{\ast }=(S_{TE}^{\ast },{\mathit{S}}_T^{\ast },{\mathit{i}}_T^{\ast },S_D^{\ast },{\mathit{i}}_D^{\ast },R_D^{\ast })\end{array}$(19) with (20) $\begin{array}{r}\{\begin{array}{l}{\displaystyle N_D^{\ast }=S_{TE}^{\ast }+{\int }_0^{\mathrm{\infty }}⟨\mathbf{1},{\mathit{i}}_T^{\ast }(a)⟩\mathrm{d}a+⟨\mathbf{1},{\mathit{S}}_T^{\ast }⟩+S_D^{\ast }}\\ {\displaystyle +{\int }_0^{\mathrm{\infty }}⟨{\mathbf{v}}_1,{\mathit{i}}_D^{\ast }(a)⟩\mathrm{d}a+{\int }_0^{\mathrm{\infty }}⟨{\mathbf{v}}_2,{\mathit{i}}_D^{\ast }(a)⟩\mathrm{d}a+R_D^{\ast }}\\ {\displaystyle 0={\theta }_D-{\mu }_D{N}_D^{\ast }-{\int }_0^{\mathrm{\infty }}{\delta }_D(a)⟨{\mathbf{v}}_1,{\mathit{i}}_D^{\ast }(a)⟩\mathrm{d}a-{\int }_0^{\mathrm{\infty }}{\delta }_C(a)⟨{\mathbf{v}}_2,{\mathit{i}}_D^{\ast }(a)⟩\mathrm{d}a}\end{array}\end{array}$(20) such that $S_{TE}^{\ast }>0$, ${\mathit{S}}_T^{\ast }\in \mathrm{int}({\mathbb{R}}_{+}^3)$, ${\mathit{i}}_T^{\ast }\ne 0_{L^1((0,+\mathrm{\infty }),{\mathbb{R}}^3)}$, $S_D^{\ast }>0$, ${\mathit{i}}_D^{\ast }\ne 0_{L^1((0,+\mathrm{\infty }),{\mathbb{R}}^2)}$, and $R_D^{\ast }>0$. By means of scaling and a suitable change of variables, we prove in Section 8 that the existence of an endemic equilibrium ${\mathcal{E}}^{\ast }$ can be reduced to the existence of a positive solution of one equation with one unknown. More precisely, we prove that there exists an endemic equilibrium ${\mathcal{E}}^{\ast }$ if and only if there exists K>0 satisfying $\begin{array}{r}\mathrm{\Delta }({\mathcal{T}}_0,K)=1\end{array}$with $\begin{array}{r}\mathrm{\Delta }({\mathcal{T}}_0,K):=\frac{1-I_D^{\diamond }({\mathcal{T}}_0,K)-(1-ϵ)R_D^{\diamond }({\mathcal{T}}_0,K)}{N_D^{\ast }(K)}f(K{\mathcal{T}}_0)\frac{1}{{\mathcal{T}}_0^{T\to D}}{\mathcal{T}}_0.\end{array}$The maps $K\to N_D^{\ast }(K)$, $K\to R_D^{\diamond }({\mathcal{T}}_0,K)$, and $K\to I_D^{\diamond }({\mathcal{T}}_0,K)$ are given by $\begin{array}{r}\{\begin{array}{l}{\displaystyle N_D^{\ast }(K):=\frac{{\theta }_D}{{\mu }_D+{\mathcal{T}}_0^{T\to D}K({\int }_0^{\mathrm{\infty }}{\delta }_D(a){\pi }_D(a)\mathrm{d}a+{\displaystyle {\int }_0^{\mathrm{\infty }}{\delta }_C(a){\pi }_C(a)\mathrm{d}a)}}}\\ {\displaystyle R_D^{\diamond }({\mathcal{T}}_0,K):={\mathcal{T}}_0^{T\to D}\frac{K{N}_D^{\ast }(K)}{ϵK{\mathcal{T}}_{0}f(K{\mathcal{T}}_0)+{\mu }_D{N}_D^{\ast }(K)}{\int }_0^{\mathrm{\infty }}\gamma (a){\pi }_D(a)\mathrm{d}a}\\ {\displaystyle I_D^{\diamond }({\mathcal{T}}_0,K)=K{\mathcal{T}}_0^{T\to D}{\int }_0^{\mathrm{\infty }}{\pi }_D(a)\mathrm{d}a+K{\mathcal{T}}_0^{T\to D}{\int }_0^{\mathrm{\infty }}{\nu }_D(a){\pi }_D(a)\mathrm{d}a{\int }_0^{\mathrm{\infty }}{\pi }_C(a)\mathrm{d}a}\end{array}\end{array}$and the function $x\to f(x)$ is defined by $\begin{array}{r}f(x):=\frac{{\theta }_D}{{\mu }_D}\sum _{k=1}^3\left(\prod _{j=1}^k\frac{{\mu }_j+{\alpha }_j}{({\mu }_j+{\alpha }_j)+x}\right)\frac{{\overline{S}}_k}{\sum _{j\in \mathcal{S}\setminus \{E\}}{\overline{S}}_j}{\mathcal{T}}_{0k}^{T\to D}\end{array}$with ${\mathcal{T}}_0^{T\to D}$ (resp. ${\mathcal{T}}_{0k}^{T\to D}$) is the vectorial (k's ticks stage vectorial) capacity given in (16(16) $\begin{array}{r}{\mathcal{T}}_0^{T\to D}=\sum _{k\in \mathcal{S}\setminus \{E\}}\underset{\begin{array}{c}\text{Proportion ticks of stage }\$k\$\end{array}}{\underbrace{\frac{{\overline{S}}_k}{\sum _{j\in \mathcal{S}\setminus \{E\}}{\overline{S}}_j}}}\times {\mathcal{T}}_{0k}^{T\to D},\end{array}$(16) ) (resp. in (17(17) $\begin{array}{r}{\mathcal{T}}_{0k}^{T\to D}=\underset{\begin{array}{c}\text{Ticks/Dogs ratio}\end{array}}{\underbrace{\frac{\sum _{j\in \mathcal{S}\setminus \{E\}}{\overline{S}}_j}{{\overline{S}}_D}}}\times \underset{\begin{array}{c}\text{Transmission from ticks at stage }k\text{ to Dogs}\end{array}}{\underbrace{{\int }_0^{\mathrm{\infty }}{\beta }_T(a)⟨\mathbf{1},e^{-\mathbf{M}a}{\mathbf{e}}_k⟩\mathrm{d}a,}}\end{array}$(17) )). Recall that in the above sum, we have used the notations $L\leftarrow 1$, $N\leftarrow 2$, and $A\leftarrow 3$ with ${\alpha }_3=0$. Furthermore, if we consider the bifurcation parameter $\begin{array}{rl}{C}_{\mathrm{bif}}^{-1}& :=\frac{1}{{\mu }_D}{\mathcal{T}}_0^{T\to D}({\int }_0^{\mathrm{\infty }}{\delta }_D(a){\pi }_D(a)\mathrm{d}a+{\int }_0^{\mathrm{\infty }}{\delta }_C(a){\pi }_C(a)\mathrm{d}a)\\ & -\sum _{k=1}^3\left(\sum _{j=1}^k\frac{1}{{\mu }_j+{\alpha }_j}\right)\frac{\frac{{\overline{S}}_k}{\sum _{j\in \mathcal{S}\setminus \{E\}}{\overline{S}}_j}{\mathcal{T}}_{0k}^{T\to D}}{{\mathcal{T}}_0^{T\to D}}\\ & -(1-ϵ)\frac{1}{{\mu }_D}{\mathcal{T}}_0^{T\to D}{\int }_0^{\mathrm{\infty }}\gamma (a){\pi }_D(a)\mathrm{d}a\\ & -{\mathcal{T}}_0^{T\to D}{\int }_0^{\mathrm{\infty }}{\pi }_D(a)\mathrm{d}a-{\mathcal{T}}_0^{T\to D}{\int }_0^{\mathrm{\infty }}{\nu }_D(a){\pi }_D(a)\mathrm{d}a{\int }_0^{\mathrm{\infty }}{\pi }_C(a)\mathrm{d}a\end{array}$then we have the following results whose proof is given in Section 8.

Theorem 3.4

Let Assumption 3.1 be satisfied. Assume in addition that the ticks' reproduction number satisfies ${\mathcal{R}}_{ST}^0>1$. Then the following properties hold true

(i)	If ${\mathcal{C}}_{\mathrm{bif}}>0$ then we have a backward bifurcation at ${\mathcal{T}}_0=1$, that is there exists an endemic equilibrium for ${\mathcal{T}}_0<1$ close to 1;
(ii)	If ${\mathcal{C}}_{\mathrm{bif}}<0$ then we have a forward bifurcation at ${\mathcal{T}}_0=1$, that is there exists an endemic equilibrium for ${\mathcal{T}}_0>1$ close to 1 and no endemic equilibrium for ${\mathcal{T}}_0<1$ close to 1.

4. Numerical simulations

In this section, we present numerical simulations to illustrate the forward and backward bifurcation results of Model (3(3) $\begin{array}{r}\{\begin{array}{l}{\displaystyle \frac{\mathrm{d}{S}_D(t)}{\mathrm{d}t}={\theta }_D-\frac{S_D(t)}{N_D(t)}{\int }_0^{\mathrm{\infty }}{\lambda }_T(t,a)\mathrm{d}a-{\mu }_D{S}_D(t),}\\ i_D(t,0)={\displaystyle \frac{S_D(t)+ϵR_D(t)}{N_D(t)}{\int }_0^{\mathrm{\infty }}{\lambda }_T(t,a)\mathrm{d}a,}\\ {\displaystyle \frac{\mathrm{\partial }{i}_D(t,a)}{\mathrm{\partial }t}+\frac{\mathrm{\partial }{i}_D(t,a)}{\mathrm{\partial }a}=-({\mu }_D+\gamma (a)+{\delta }_D(a)+{\nu }_D(a))i_D(t,a),}\\ {\displaystyle \frac{\mathrm{\partial }{c}_D(t,a)}{\mathrm{\partial }t}+\frac{\mathrm{\partial }{c}_D(t,a)}{\mathrm{\partial }a}=-({\mu }_D+{\delta }_C(a))c_D(t,a)}\\ c_D(t,0)={\displaystyle {\int }_0^{\mathrm{\infty }}{\nu }_D(a)i_D(t,a)\mathrm{d}a,}\\ {\displaystyle \frac{\mathrm{d}{R}_D(t)}{\mathrm{d}t}={\int }_0^{\mathrm{\infty }}\gamma (a)i_D(t,a)\mathrm{d}a-ϵ\frac{R_D(t)}{N_D(t)}{\int }_0^{\mathrm{\infty }}{\lambda }_T(t,a)\mathrm{d}a-{\mu }_D{R}_D(t).}\end{array}\end{array}$(3) )–(4(4) $\begin{array}{r}\{\begin{array}{l}{\displaystyle \frac{\mathrm{d}{S}_{TE}(t)}{\mathrm{d}t}=r_E(1-\frac{S_{TE}(t)}{K})N_{TA}(t)-({\alpha }_E+{\mu }_E)S_{TE}(t),}\\ {\displaystyle \frac{\mathrm{d}{S}_{TL}(t)}{\mathrm{d}t}={\alpha }_E{S}_{TE}(t)-\frac{S_{TL}(t)}{N_D(t)}{\lambda }_D(t)-({\alpha }_L+{\mu }_L)S_{TL}(t),}\\ {\displaystyle \frac{\mathrm{d}{S}_{TN}(t)}{\mathrm{d}t}={\alpha }_L{S}_{TL}(t)-\frac{S_{TN}(t)}{N_D(t)}{\lambda }_D(t)-({\alpha }_N+{\mu }_N)S_{TN}(t),}\\ {\displaystyle \frac{\mathrm{d}{S}_{TA}(t)}{\mathrm{d}t}={\alpha }_N{S}_{TN}(t)-\frac{S_{TA}(t)}{N_D(t)}{\lambda }_D(t)-{\mu }_A{S}_{TA}(t),}\\ \text{and for }k\in \mathcal{S}\setminus \{E\},\\ i_{Tk}(t,0)={\displaystyle \frac{S_{Tk}(t)}{N_D(t)}{\lambda }_D(t),}\\ {\displaystyle \frac{\mathrm{\partial }{i}_{Tk}(t,a)}{\mathrm{\partial }t}+\frac{\mathrm{\partial }{i}_{Tk}(t,a)}{\mathrm{\partial }a}=-({\mu }_k+{\alpha }_k)i_{Tk}(t,a)}\end{array}\end{array}$(4) ). These simulations were conducted using finite volume numerical schemes implemented with the R software (http://www.r-project.org/). The values of the parameters used are given in Table 1. In all the figures, we have used the initial conditions $S_D(0)=10$, $i_D(0,a)=0.03amax(0,100-a)$, $c_D(0,a)=0.01amax(0,100-a)$, $R_D(0)=0$, $S_{TE}(0)=1$, $S_{TL}(0)=10$, $S_{TN}(0)=0$, $S_{TA}(0)=2$, $i_{TL}(0,a)=0.01$, $i_{TN}(0,a)=0.3$, $i_{TA}(0,a)=1$. Given the parameter values in Table 1, the ticks' reproduction number ${\mathcal{R}}_{ST}^0\approx 66.193$.

A forward bifurcation occurs at ${\mathcal{T}}_0=1$ (Theorem 3.4), which means that whenever ${\mathcal{T}}_0<1$, then the disease-free equilibrium ${\mathcal{E}}_0$ is locally asymptotically stable (Theorem 3.3) and no endemic equilibrium exists. Asymptotically, the disease go extinct (Figure 2(b), where ${\mathcal{T}}_0\approx 0.992$). However, if ${\mathcal{T}}_0>1$, then ${\mathcal{E}}_0$ is unstable (Theorem 3.3) and an endemic equilibrium exists if $C_{\mathrm{bif}}<0$ (Theorem 3.4). The disease is asymptotically persistent and the solution converges to the endemic equilibrium (Figure 2(a), where ${\mathcal{T}}_0\approx 1.103$ and $C_{\mathrm{bif}}\approx -9.4\times {10}^{-4}$).

Figure 2. Forward bifurcation diagram with convergence either to (a) the endemic equilibrium when ${\mathcal{T}}_0\approx 1.103>1$; (b) the disease-free equilibrium when ${\mathcal{T}}_0\approx 0.992<1$. Here $r_E=1$, ${\delta }_D(a)=1.735\times {10}^{-4}$, ${\delta }_C(a)=3.47\times {10}^{-4}$, ${\beta }_D^2(a)=0.038$, ${\beta }_T(a)=0.15504$, and $C_{\mathrm{bif}}=-9.4\times {10}^{-4}$ are fixed. The probability of infection ${\beta }_D^1(a)$ is considered constant with: (a) ${\beta }_D^1(a)=0.076$ and (b) ${\beta }_D^1(a)=0.0506$. The other parameters are given by Table 1.

A backward bifurcation occurs at ${\mathcal{T}}_0=1$. This means that whenever ${\mathcal{T}}_0>1$, the disease-free equilibrium is unstable and there exists a unique endemic equilibrium. In such a situation, the solutions converge asymptotically to this endemic equilibrium (Figures 3(a)), where ${\mathcal{T}}_0\approx 1.318$). By contrast to the forward bifurcation (Figures 2 and 4(a)), when ${\mathcal{T}}_0<1$ and $C_{\mathrm{bif}}\approx 5.21\times {10}^{-4}>0$ there exists an endemic equilibrium (Theorem 3.4). Furthermore, there exists a threshold ${\mathcal{T}}_0^{\ast }\approx 0.422<1$ such that (i) for ${\mathcal{T}}_0<{\mathcal{T}}_0^{\ast }$, there is no endemic equilibrium and the solutions converge to the disease-free equilibrium (Figures 3(b) and 4(b), where ${\mathcal{T}}_0\approx 0.415<{\mathcal{T}}_0^{\ast }$) and (ii) for ${\mathcal{T}}_0^{\ast }<{\mathcal{T}}_0<1$, there exist an endemic equilibrium such that, depending on the initial condition, the solution can converge to this endemic equilibrium (Figures 3(c) and 4(b), where ${\mathcal{T}}_0^{\ast }<{\mathcal{T}}_0\approx 0.843$).

Figure 3. Backward bifurcation diagram with convergence either to (a) the endemic equilibrium when ${\mathcal{T}}_0\approx 1.318>1$; (b) the disease-free equilibrium when ${\mathcal{T}}_0\approx 0.415<{\mathcal{T}}_0^{\ast }<1$; (c) the endemic equilibrium when ${\mathcal{T}}_0^{\ast }<{\mathcal{T}}_0\approx 0.843<1$. Here $r_E=1$, ${\delta }_D(a)=0.1735$, ${\delta }_C(a)=0.347$, ${\beta }_D^2(a)=0.608$, ${\beta }_T(a)=0.15504$, and $C_{\mathrm{bif}}=5.21\times {10}^{-4}$ are fixed. The probability of infection ${\beta }_D^1(a)$ is considered constant with: (a) ${\beta }_D^1(a)=0.76$, (b) ${\beta }_D^1(a)=0.304$, and (c) ${\beta }_D^1(a)=0.456$. The other parameters are given by Table 1.

Figure 4. Forward and backward bifurcation diagram. (a) Forward bifurcation diagram with $C_{\mathrm{bif}}=-9.4\times {10}^{-4}$ where all the parameters are the same as in Figure 2. (b) Backward bifurcation diagram with $C_{\mathrm{bif}}=5.21\times {10}^{-4}$ and the parameters are the same as in Figure 3. Note that $C_{\mathrm{bif}}$ does not depend on the parameter ${\beta }_D^1$.

5. Discussion and conclusions

5.1. Discussion

In this work, we have developed a novel partial differential equation (PDE) model that extends classical epidemiological models proposed for tick-borne diseases (see [29,30]). The model proposed here incorporates the ability to accurately account for the different ticks' developmental stages (discrete variable) as well as the variations in infectiousness over the course of infection (continuous variable). In our study, we established the mathematical well-posedness of the model using the integrated semigroups theory. However, the presence of singularity in the force of infection whenever the total dogs' population is zero introduces a complication to the analysis. Note that this is a mathematical construct, since without dogs the model would not hold nor make sense.

We derived an explicit formula for the reproduction number ${\mathcal{T}}_0$, extending the classical formula. We identified two possible behaviours around ${\mathcal{T}}_0=1$. The first scenario involves a forward bifurcation, indicating that an epidemic can only occur if ${\mathcal{T}}_0>1$. In the second scenario, a backward bifurcation is observed, where an epidemic can arise if ${\mathcal{T}}_0<1$, provided that ${\mathcal{T}}_0$ is sufficiently close to 1. These findings have significant epidemiological implications, especially in an endemic region where controlling the epidemic is of importance. In the forward bifurcation case, simply reducing the reproduction number ${\mathcal{T}}_0$ below 1 is sufficient to halt the epidemic. However, in the backward bifurcation scenario, the number must be reduced below a second threshold, denoted as ${\mathcal{T}}_0^{\ast }$.

The epidemiological implication of the backward bifurcation phenomenon is that the classical requirement of having the basic reproduction number (${\mathcal{R}}_0$) less than one is no longer sufficient to ensure effective disease eradication or elimination [38]. This implies that disease can invade to a relatively high endemic level once the reproduction number ${\mathcal{R}}_0$ is more than one; decreasing ${\mathcal{R}}_0$ below one do not necessarily make the disease disappear, see Figure 3(b) and Figure 2 in [38]. Backward bifurcation has been shown to occur in several vector-borne disease models [37,39–45]. For instance Garba et al. [44] showed for a dengue model the possibility of backward bifurcation where a locally stable disease-free equilibrium coexists with a locally stable endemic equilibrium. Also, the age-structure malaria and Chikungunya models in [37,39,42,43] were shown to exhibit backward bifurcation when ${\mathcal{R}}_0$ is less than one. This phenomenon is equally possible in the absence of any age-structure [39,42,43]. Backward bifurcation in these vector-borne disease models is induced by disease induced mortality. Furthermore, backward bifurcation can equally be induced by several other factors like vaccination [46–49], exogenous re-infection which are known to occur in tuberculosis models [50–52], and cross-immunity for instance in a two strain influenza transmission model [53]. Other epidemiological mechanism like differential susceptibility in risk-structured models could also cause backward bifurcation in disease transmission models [38]. We should note that the age-structured Ehrlichia chaffeensis model (3(3) $\begin{array}{r}\{\begin{array}{l}{\displaystyle \frac{\mathrm{d}{S}_D(t)}{\mathrm{d}t}={\theta }_D-\frac{S_D(t)}{N_D(t)}{\int }_0^{\mathrm{\infty }}{\lambda }_T(t,a)\mathrm{d}a-{\mu }_D{S}_D(t),}\\ i_D(t,0)={\displaystyle \frac{S_D(t)+ϵR_D(t)}{N_D(t)}{\int }_0^{\mathrm{\infty }}{\lambda }_T(t,a)\mathrm{d}a,}\\ {\displaystyle \frac{\mathrm{\partial }{i}_D(t,a)}{\mathrm{\partial }t}+\frac{\mathrm{\partial }{i}_D(t,a)}{\mathrm{\partial }a}=-({\mu }_D+\gamma (a)+{\delta }_D(a)+{\nu }_D(a))i_D(t,a),}\\ {\displaystyle \frac{\mathrm{\partial }{c}_D(t,a)}{\mathrm{\partial }t}+\frac{\mathrm{\partial }{c}_D(t,a)}{\mathrm{\partial }a}=-({\mu }_D+{\delta }_C(a))c_D(t,a)}\\ c_D(t,0)={\displaystyle {\int }_0^{\mathrm{\infty }}{\nu }_D(a)i_D(t,a)\mathrm{d}a,}\\ {\displaystyle \frac{\mathrm{d}{R}_D(t)}{\mathrm{d}t}={\int }_0^{\mathrm{\infty }}\gamma (a)i_D(t,a)\mathrm{d}a-ϵ\frac{R_D(t)}{N_D(t)}{\int }_0^{\mathrm{\infty }}{\lambda }_T(t,a)\mathrm{d}a-{\mu }_D{R}_D(t).}\end{array}\end{array}$(3) )–(4(4) $\begin{array}{r}\{\begin{array}{l}{\displaystyle \frac{\mathrm{d}{S}_{TE}(t)}{\mathrm{d}t}=r_E(1-\frac{S_{TE}(t)}{K})N_{TA}(t)-({\alpha }_E+{\mu }_E)S_{TE}(t),}\\ {\displaystyle \frac{\mathrm{d}{S}_{TL}(t)}{\mathrm{d}t}={\alpha }_E{S}_{TE}(t)-\frac{S_{TL}(t)}{N_D(t)}{\lambda }_D(t)-({\alpha }_L+{\mu }_L)S_{TL}(t),}\\ {\displaystyle \frac{\mathrm{d}{S}_{TN}(t)}{\mathrm{d}t}={\alpha }_L{S}_{TL}(t)-\frac{S_{TN}(t)}{N_D(t)}{\lambda }_D(t)-({\alpha }_N+{\mu }_N)S_{TN}(t),}\\ {\displaystyle \frac{\mathrm{d}{S}_{TA}(t)}{\mathrm{d}t}={\alpha }_N{S}_{TN}(t)-\frac{S_{TA}(t)}{N_D(t)}{\lambda }_D(t)-{\mu }_A{S}_{TA}(t),}\\ \text{and for }k\in \mathcal{S}\setminus \{E\},\\ i_{Tk}(t,0)={\displaystyle \frac{S_{Tk}(t)}{N_D(t)}{\lambda }_D(t),}\\ {\displaystyle \frac{\mathrm{\partial }{i}_{Tk}(t,a)}{\mathrm{\partial }t}+\frac{\mathrm{\partial }{i}_{Tk}(t,a)}{\mathrm{\partial }a}=-({\mu }_k+{\alpha }_k)i_{Tk}(t,a)}\end{array}\end{array}$(4) ) include disease induce death (${\delta }_D$ and ${\delta }_C$) in the dog population and hence the source of the backward bifurcation.

5.2. Conclusion

In conclusion, we have developed a novel partial differential equation model of Ehrlichia chaffeensis transmission dynamics in dogs. The model incorporates the different developmental life stages of ticks (discrete variable) as well as the duration of infection (continuous variable). The following results were obtained from our theoretical analysis and numerical simulations:

1. The developed model is well-posed;

2. The model always exhibits a disease-free equilibrium along with an endemic equilibrium;

3. The model has a reproduction number, denoted as ${\mathcal{T}}_0$;

4. A necessary and sufficient condition for the bifurcation of an endemic equilibrium was established using semigroup approach;

5. A bifurcation (forward or backward), can occur at ${\mathcal{T}}_0=1$.

6. The disease-free equilibrium of (3)–(4)

The dynamic of ticks in the absence of infection is governed by the following system of equations (21) $\begin{array}{r}\{\begin{array}{l}{\displaystyle \frac{\mathrm{d}{S}_{TE}(t)}{\mathrm{d}t}=r_E(1-\frac{S_{TE}(t)}{K})S_{TA}(t)-({\mu }_E+{\alpha }_E)S_{TE}(t)}\\ {\displaystyle \frac{\mathrm{d}{S}_{TL}(t)}{\mathrm{d}t}={\alpha }_E{S}_{TE}(t)-({\mu }_L+{\alpha }_L)S_{TL}(t)}\\ {\displaystyle \frac{\mathrm{d}{S}_{TN}(t)}{\mathrm{d}t}={\alpha }_L{S}_{TL}(t)-({\mu }_N+{\alpha }_N)S_{TN}(t)}\\ {\displaystyle \frac{\mathrm{d}{S}_{TA}(t)}{\mathrm{d}t}={\alpha }_N{S}_{TN}(t)-{\mu }_A{S}_{TA}(t)}\end{array}\end{array}$(21) with initial condition $S_{TE}(0)=S_{TE0}\ge 0$, $S_{TL}(0)=S_{TL0}\ge 0$, $S_{TN}(0)=S_{TN0}\ge 0$, and $S_{TA}(0)=S_{TA0}\ge 0$. It is easy to see that $\{0\}\times {\mathbb{R}}_{+}^3$ is invariant with respect to the ordinary differential equation so that if $S_{TE0}=0$ then the larval, nymphal, and adult stages exponentially goes to 0. From this, we have that a necessary condition for the ticks to persist is $S_{TE0}>0$. Moreover, by straightforward computations, we find that (21(21) $\begin{array}{r}\{\begin{array}{l}{\displaystyle \frac{\mathrm{d}{S}_{TE}(t)}{\mathrm{d}t}=r_E(1-\frac{S_{TE}(t)}{K})S_{TA}(t)-({\mu }_E+{\alpha }_E)S_{TE}(t)}\\ {\displaystyle \frac{\mathrm{d}{S}_{TL}(t)}{\mathrm{d}t}={\alpha }_E{S}_{TE}(t)-({\mu }_L+{\alpha }_L)S_{TL}(t)}\\ {\displaystyle \frac{\mathrm{d}{S}_{TN}(t)}{\mathrm{d}t}={\alpha }_L{S}_{TL}(t)-({\mu }_N+{\alpha }_N)S_{TN}(t)}\\ {\displaystyle \frac{\mathrm{d}{S}_{TA}(t)}{\mathrm{d}t}={\alpha }_N{S}_{TN}(t)-{\mu }_A{S}_{TA}(t)}\end{array}\end{array}$(21) ) has a unique positive stationary solution if and only if the ticks reproduction number ${\mathcal{R}}_{ST}^0$ defined by (12(12) $\begin{array}{r}{\mathcal{R}}_{ST}^0=\frac{r_E}{{\mu }_A}\prod _{k\in \mathcal{S}}\frac{{\alpha }_k}{{\alpha }_k+{\mu }_k}=r_E\frac{{\alpha }_E}{{\alpha }_E+{\mu }_E}\frac{{\alpha }_L}{{\alpha }_L+{\mu }_L}\frac{{\alpha }_N}{{\alpha }_N+{\mu }_N}\frac{1}{{\mu }_A}.\end{array}$(12) ) satisfies ${\mathcal{R}}_{ST}^0>1$. By straightforward computations, we have that the equilibrium of system (21(21) $\begin{array}{r}\{\begin{array}{l}{\displaystyle \frac{\mathrm{d}{S}_{TE}(t)}{\mathrm{d}t}=r_E(1-\frac{S_{TE}(t)}{K})S_{TA}(t)-({\mu }_E+{\alpha }_E)S_{TE}(t)}\\ {\displaystyle \frac{\mathrm{d}{S}_{TL}(t)}{\mathrm{d}t}={\alpha }_E{S}_{TE}(t)-({\mu }_L+{\alpha }_L)S_{TL}(t)}\\ {\displaystyle \frac{\mathrm{d}{S}_{TN}(t)}{\mathrm{d}t}={\alpha }_L{S}_{TL}(t)-({\mu }_N+{\alpha }_N)S_{TN}(t)}\\ {\displaystyle \frac{\mathrm{d}{S}_{TA}(t)}{\mathrm{d}t}={\alpha }_N{S}_{TN}(t)-{\mu }_A{S}_{TA}(t)}\end{array}\end{array}$(21) ) is given by (22) $\begin{array}{r}\{\begin{array}{l}{\displaystyle {\overline{S}}_{TE}=K\frac{{\mathcal{R}}_{ST}^0-1}{{\mathcal{R}}_{ST}^0},}\\ {\displaystyle {\overline{S}}_{TL}=\frac{{\alpha }_E}{{\mu }_L+{\alpha }_L}{\overline{S}}_{TE},}\\ {\displaystyle {\overline{S}}_{TN}=\frac{{\alpha }_L}{{\mu }_N+{\alpha }_N}\frac{{\alpha }_E}{{\mu }_L+{\alpha }_L}{\overline{S}}_{TE},}\\ {\displaystyle {\overline{S}}_{TA}=\frac{{\alpha }_N}{{\mu }_A}\frac{{\alpha }_L}{{\mu }_N+{\alpha }_N}\frac{{\alpha }_E}{{\mu }_L+{\alpha }_L}{\overline{S}}_{TE}.}\end{array}\end{array}$(22) Note that in accordance with the compact formulation (7(7) $\begin{array}{r}\{\begin{array}{l}{\displaystyle \frac{\mathrm{d}{S}_D(t)}{\mathrm{d}t}={\theta }_D-\frac{S_D(t)}{N_D(t)}{\int }_0^{\mathrm{\infty }}{\beta }_T(a)⟨\mathbf{1},{\mathit{i}}_T(t,a)⟩\mathrm{d}a-{\mu }_D{S}_D(t),}\\ {\displaystyle \frac{\mathrm{\partial }{\mathit{i}}_D(t,a)}{\mathrm{\partial }t}+\frac{\mathrm{\partial }{\mathit{i}}_D(t,a)}{\mathrm{\partial }a}=-({\mu }_D+\mathit{\theta }(a)){\mathit{i}}_D(t,a),}\\ {\mathit{i}}_D(t,0)={\displaystyle \frac{S_D(t)+ϵR_D(t)}{N_D(t)}{\int }_0^{\mathrm{\infty }}{\beta }_T(a)⟨\mathbf{1},{\mathit{i}}_T(t,a)⟩\mathrm{d}a{\mathbf{v}}_1}\\ {\displaystyle +{\int }_0^{\mathrm{\infty }}{\nu }_D(a)⟨{\mathbf{v}}_1,{\mathit{i}}_D(t,a)⟩\mathrm{d}a{\mathbf{v}}_2,}\\ {\displaystyle \frac{\mathrm{d}{R}_D(t)}{\mathrm{d}t}={\int }_0^{\mathrm{\infty }}\gamma (a)⟨{\mathbf{v}}_1,{\mathit{i}}_D(t,a)⟩\mathrm{d}a-ϵ\frac{R_D(t)}{N_D(t)}{\int }_0^{\mathrm{\infty }}{\beta }_T(a)⟨\mathbf{1},{\mathit{i}}_T(t,a)⟩\mathrm{d}a-{\mu }_D{R}_D(t),}\end{array}\end{array}$(7) )–(8(8) $\begin{array}{r}\{\begin{array}{l}{\displaystyle \frac{\mathrm{d}{S}_{TE}(t)}{\mathrm{d}t}=r_E(1-\frac{S_{TE}(t)}{K}){\displaystyle (⟨{\mathit{S}}_T(t),{\mathbf{e}}_3⟩+{\int }_0^{\mathrm{\infty }}⟨{\mathit{i}}_T(t,a),{\mathbf{e}}_3⟩\mathrm{d}a)-({\alpha }_E+{\mu }_E)S_{TE}(t),}}\\ {\displaystyle \frac{\mathrm{d}{\mathit{S}}_T(t)}{\mathrm{d}t}={\alpha }_E{S}_{TE}(t){\mathbf{e}}_1-\mathbf{M}{\mathit{S}}_T(t)-\frac{{\mathit{S}}_T(t)}{N_D(t)}{\int }_0^{\mathrm{\infty }}⟨{\mathit{\beta }}_D(a),{\mathit{i}}_D(t,a)⟩\mathrm{d}a,}\\ {\displaystyle \frac{\mathrm{\partial }{\mathit{i}}_T(t,a)}{\mathrm{\partial }t}+\frac{\mathrm{\partial }{\mathit{i}}_T(t,a)}{\mathrm{\partial }a}=-\mathbf{M}{\mathit{i}}_T(t,a),}\\ {\mathit{i}}_T(t,0)={\displaystyle \frac{{\mathit{S}}_T(t)}{N_D(t)}{\int }_0^{\mathrm{\infty }}⟨{\mathit{\beta }}_D(a),{\mathit{i}}_D(t,a)⟩\mathrm{d}a,}\end{array}\end{array}$(8) ), the above expressions of ${\mathcal{R}}_{ST}^0$ and ${\overline{S}}_{Tk}$, with $k\in \mathcal{S}$, can be rewritten as (23) $\begin{array}{r}{\mathcal{R}}_{ST}^0:=\frac{r_E{\alpha }_E}{{\alpha }_E+{\mu }_E}⟨{\mathbf{M}}^{-1}{\mathbf{e}}_1,{\mathbf{e}}_3⟩,\end{array}$(23) and (24) $\begin{array}{r}\{\begin{array}{l}{\overline{S}}_{TE}={\displaystyle K(1-\frac{1}{{\mathcal{R}}_{ST}^0}),}\\ {\displaystyle {\overline{\mathit{S}}}_T={\alpha }_E{\overline{S}}_{TE}\sum _{k=1}^3{\mathbf{e}}_k\prod _{j=1}^k\frac{{\alpha }_{j-1}}{{\mu }_j+{\alpha }_j},\text{with }{\alpha }_0=1,{\alpha }_3=0,}\end{array}\end{array}$(24) and using the correspondence $L\leftarrow 1$, $N\leftarrow 2$, and $A\leftarrow 3$.

7. Threshold and proof of Theorem 3.3

In order to obtain the threshold of (7(7) $\begin{array}{r}\{\begin{array}{l}{\displaystyle \frac{\mathrm{d}{S}_D(t)}{\mathrm{d}t}={\theta }_D-\frac{S_D(t)}{N_D(t)}{\int }_0^{\mathrm{\infty }}{\beta }_T(a)⟨\mathbf{1},{\mathit{i}}_T(t,a)⟩\mathrm{d}a-{\mu }_D{S}_D(t),}\\ {\displaystyle \frac{\mathrm{\partial }{\mathit{i}}_D(t,a)}{\mathrm{\partial }t}+\frac{\mathrm{\partial }{\mathit{i}}_D(t,a)}{\mathrm{\partial }a}=-({\mu }_D+\mathit{\theta }(a)){\mathit{i}}_D(t,a),}\\ {\mathit{i}}_D(t,0)={\displaystyle \frac{S_D(t)+ϵR_D(t)}{N_D(t)}{\int }_0^{\mathrm{\infty }}{\beta }_T(a)⟨\mathbf{1},{\mathit{i}}_T(t,a)⟩\mathrm{d}a{\mathbf{v}}_1}\\ {\displaystyle +{\int }_0^{\mathrm{\infty }}{\nu }_D(a)⟨{\mathbf{v}}_1,{\mathit{i}}_D(t,a)⟩\mathrm{d}a{\mathbf{v}}_2,}\\ {\displaystyle \frac{\mathrm{d}{R}_D(t)}{\mathrm{d}t}={\int }_0^{\mathrm{\infty }}\gamma (a)⟨{\mathbf{v}}_1,{\mathit{i}}_D(t,a)⟩\mathrm{d}a-ϵ\frac{R_D(t)}{N_D(t)}{\int }_0^{\mathrm{\infty }}{\beta }_T(a)⟨\mathbf{1},{\mathit{i}}_T(t,a)⟩\mathrm{d}a-{\mu }_D{R}_D(t),}\end{array}\end{array}$(7) )–(8(8) $\begin{array}{r}\{\begin{array}{l}{\displaystyle \frac{\mathrm{d}{S}_{TE}(t)}{\mathrm{d}t}=r_E(1-\frac{S_{TE}(t)}{K}){\displaystyle (⟨{\mathit{S}}_T(t),{\mathbf{e}}_3⟩+{\int }_0^{\mathrm{\infty }}⟨{\mathit{i}}_T(t,a),{\mathbf{e}}_3⟩\mathrm{d}a)-({\alpha }_E+{\mu }_E)S_{TE}(t),}}\\ {\displaystyle \frac{\mathrm{d}{\mathit{S}}_T(t)}{\mathrm{d}t}={\alpha }_E{S}_{TE}(t){\mathbf{e}}_1-\mathbf{M}{\mathit{S}}_T(t)-\frac{{\mathit{S}}_T(t)}{N_D(t)}{\int }_0^{\mathrm{\infty }}⟨{\mathit{\beta }}_D(a),{\mathit{i}}_D(t,a)⟩\mathrm{d}a,}\\ {\displaystyle \frac{\mathrm{\partial }{\mathit{i}}_T(t,a)}{\mathrm{\partial }t}+\frac{\mathrm{\partial }{\mathit{i}}_T(t,a)}{\mathrm{\partial }a}=-\mathbf{M}{\mathit{i}}_T(t,a),}\\ {\mathit{i}}_T(t,0)={\displaystyle \frac{{\mathit{S}}_T(t)}{N_D(t)}{\int }_0^{\mathrm{\infty }}⟨{\mathit{\beta }}_D(a),{\mathit{i}}_D(t,a)⟩\mathrm{d}a,}\end{array}\end{array}$(8) ) that determine disease invasion dynamics in a completely susceptible population of ticks and dogs, we consider the linearized equation to (7(7) $\begin{array}{r}\{\begin{array}{l}{\displaystyle \frac{\mathrm{d}{S}_D(t)}{\mathrm{d}t}={\theta }_D-\frac{S_D(t)}{N_D(t)}{\int }_0^{\mathrm{\infty }}{\beta }_T(a)⟨\mathbf{1},{\mathit{i}}_T(t,a)⟩\mathrm{d}a-{\mu }_D{S}_D(t),}\\ {\displaystyle \frac{\mathrm{\partial }{\mathit{i}}_D(t,a)}{\mathrm{\partial }t}+\frac{\mathrm{\partial }{\mathit{i}}_D(t,a)}{\mathrm{\partial }a}=-({\mu }_D+\mathit{\theta }(a)){\mathit{i}}_D(t,a),}\\ {\mathit{i}}_D(t,0)={\displaystyle \frac{S_D(t)+ϵR_D(t)}{N_D(t)}{\int }_0^{\mathrm{\infty }}{\beta }_T(a)⟨\mathbf{1},{\mathit{i}}_T(t,a)⟩\mathrm{d}a{\mathbf{v}}_1}\\ {\displaystyle +{\int }_0^{\mathrm{\infty }}{\nu }_D(a)⟨{\mathbf{v}}_1,{\mathit{i}}_D(t,a)⟩\mathrm{d}a{\mathbf{v}}_2,}\\ {\displaystyle \frac{\mathrm{d}{R}_D(t)}{\mathrm{d}t}={\int }_0^{\mathrm{\infty }}\gamma (a)⟨{\mathbf{v}}_1,{\mathit{i}}_D(t,a)⟩\mathrm{d}a-ϵ\frac{R_D(t)}{N_D(t)}{\int }_0^{\mathrm{\infty }}{\beta }_T(a)⟨\mathbf{1},{\mathit{i}}_T(t,a)⟩\mathrm{d}a-{\mu }_D{R}_D(t),}\end{array}\end{array}$(7) )–(8(8) $\begin{array}{r}\{\begin{array}{l}{\displaystyle \frac{\mathrm{d}{S}_{TE}(t)}{\mathrm{d}t}=r_E(1-\frac{S_{TE}(t)}{K}){\displaystyle (⟨{\mathit{S}}_T(t),{\mathbf{e}}_3⟩+{\int }_0^{\mathrm{\infty }}⟨{\mathit{i}}_T(t,a),{\mathbf{e}}_3⟩\mathrm{d}a)-({\alpha }_E+{\mu }_E)S_{TE}(t),}}\\ {\displaystyle \frac{\mathrm{d}{\mathit{S}}_T(t)}{\mathrm{d}t}={\alpha }_E{S}_{TE}(t){\mathbf{e}}_1-\mathbf{M}{\mathit{S}}_T(t)-\frac{{\mathit{S}}_T(t)}{N_D(t)}{\int }_0^{\mathrm{\infty }}⟨{\mathit{\beta }}_D(a),{\mathit{i}}_D(t,a)⟩\mathrm{d}a,}\\ {\displaystyle \frac{\mathrm{\partial }{\mathit{i}}_T(t,a)}{\mathrm{\partial }t}+\frac{\mathrm{\partial }{\mathit{i}}_T(t,a)}{\mathrm{\partial }a}=-\mathbf{M}{\mathit{i}}_T(t,a),}\\ {\mathit{i}}_T(t,0)={\displaystyle \frac{{\mathit{S}}_T(t)}{N_D(t)}{\int }_0^{\mathrm{\infty }}⟨{\mathit{\beta }}_D(a),{\mathit{i}}_D(t,a)⟩\mathrm{d}a,}\end{array}\end{array}$(8) ) at the disease-free equilibrium ${\mathcal{E}}_0$. More precisely, we linearize the infective compartments of (7(7) $\begin{array}{r}\{\begin{array}{l}{\displaystyle \frac{\mathrm{d}{S}_D(t)}{\mathrm{d}t}={\theta }_D-\frac{S_D(t)}{N_D(t)}{\int }_0^{\mathrm{\infty }}{\beta }_T(a)⟨\mathbf{1},{\mathit{i}}_T(t,a)⟩\mathrm{d}a-{\mu }_D{S}_D(t),}\\ {\displaystyle \frac{\mathrm{\partial }{\mathit{i}}_D(t,a)}{\mathrm{\partial }t}+\frac{\mathrm{\partial }{\mathit{i}}_D(t,a)}{\mathrm{\partial }a}=-({\mu }_D+\mathit{\theta }(a)){\mathit{i}}_D(t,a),}\\ {\mathit{i}}_D(t,0)={\displaystyle \frac{S_D(t)+ϵR_D(t)}{N_D(t)}{\int }_0^{\mathrm{\infty }}{\beta }_T(a)⟨\mathbf{1},{\mathit{i}}_T(t,a)⟩\mathrm{d}a{\mathbf{v}}_1}\\ {\displaystyle +{\int }_0^{\mathrm{\infty }}{\nu }_D(a)⟨{\mathbf{v}}_1,{\mathit{i}}_D(t,a)⟩\mathrm{d}a{\mathbf{v}}_2,}\\ {\displaystyle \frac{\mathrm{d}{R}_D(t)}{\mathrm{d}t}={\int }_0^{\mathrm{\infty }}\gamma (a)⟨{\mathbf{v}}_1,{\mathit{i}}_D(t,a)⟩\mathrm{d}a-ϵ\frac{R_D(t)}{N_D(t)}{\int }_0^{\mathrm{\infty }}{\beta }_T(a)⟨\mathbf{1},{\mathit{i}}_T(t,a)⟩\mathrm{d}a-{\mu }_D{R}_D(t),}\end{array}\end{array}$(7) )–(8(8) $\begin{array}{r}\{\begin{array}{l}{\displaystyle \frac{\mathrm{d}{S}_{TE}(t)}{\mathrm{d}t}=r_E(1-\frac{S_{TE}(t)}{K}){\displaystyle (⟨{\mathit{S}}_T(t),{\mathbf{e}}_3⟩+{\int }_0^{\mathrm{\infty }}⟨{\mathit{i}}_T(t,a),{\mathbf{e}}_3⟩\mathrm{d}a)-({\alpha }_E+{\mu }_E)S_{TE}(t),}}\\ {\displaystyle \frac{\mathrm{d}{\mathit{S}}_T(t)}{\mathrm{d}t}={\alpha }_E{S}_{TE}(t){\mathbf{e}}_1-\mathbf{M}{\mathit{S}}_T(t)-\frac{{\mathit{S}}_T(t)}{N_D(t)}{\int }_0^{\mathrm{\infty }}⟨{\mathit{\beta }}_D(a),{\mathit{i}}_D(t,a)⟩\mathrm{d}a,}\\ {\displaystyle \frac{\mathrm{\partial }{\mathit{i}}_T(t,a)}{\mathrm{\partial }t}+\frac{\mathrm{\partial }{\mathit{i}}_T(t,a)}{\mathrm{\partial }a}=-\mathbf{M}{\mathit{i}}_T(t,a),}\\ {\mathit{i}}_T(t,0)={\displaystyle \frac{{\mathit{S}}_T(t)}{N_D(t)}{\int }_0^{\mathrm{\infty }}⟨{\mathit{\beta }}_D(a),{\mathit{i}}_D(t,a)⟩\mathrm{d}a,}\end{array}\end{array}$(8) ) around ${\mathcal{E}}_0$ to obtain the following system (25) $\begin{array}{r}\{\begin{array}{l}{\displaystyle \frac{\mathrm{\partial }{\mathit{i}}_D(t,a)}{\mathrm{\partial }t}+\frac{\mathrm{\partial }{\mathit{i}}_D(t,a)}{\mathrm{\partial }a}=-({\mu }_D+\mathit{\theta }(a)){\mathit{i}}_D(t,a)}\\ {\mathit{i}}_D(t,0)={\displaystyle {\int }_0^{\mathrm{\infty }}{\beta }_T(a)⟨\mathbf{1},{\mathit{i}}_T(t,a)⟩\mathrm{d}a{\mathbf{v}}_1}\\ {\displaystyle +{\int }_0^{\mathrm{\infty }}{\nu }_D(a)⟨{\mathbf{v}}_1,{\mathit{i}}_D(t,a)⟩\mathrm{d}a{\mathbf{v}}_2}\end{array}\end{array}$(25) and (26) $\begin{array}{r}\{\begin{array}{l}\frac{\mathrm{\partial }{\mathit{i}}_T(t,a)}{\mathrm{\partial }t}+\frac{\mathrm{\partial }{\mathit{i}}_T(t,a)}{\mathrm{\partial }a}=-\mathbf{M}{\mathit{i}}_T(t,a)\\ {\mathit{i}}_T(t,0)={\displaystyle \frac{{\overline{\mathit{S}}}_T}{{\overline{S}}_D}{\int }_0^{\mathrm{\infty }}⟨{\mathit{\beta }}_D(a),{\mathit{i}}_D(t,a)⟩\mathrm{d}a}\end{array}\end{array}$(26) whose initial conditions satisfy ${\mathit{i}}_D(0,\cdot )\in L^1((0,+\mathrm{\infty }),{\mathbb{R}}^2)$ and ${\mathit{i}}_D(0,\cdot )\in L^1((0,+\mathrm{\infty }),{\mathbb{R}}^3)$. In order to define the threshold and study the local asymptotic stability of the disease-free equilibrium, we will make use of the semigroup approach by reformulating (25(25) $\begin{array}{r}\{\begin{array}{l}{\displaystyle \frac{\mathrm{\partial }{\mathit{i}}_D(t,a)}{\mathrm{\partial }t}+\frac{\mathrm{\partial }{\mathit{i}}_D(t,a)}{\mathrm{\partial }a}=-({\mu }_D+\mathit{\theta }(a)){\mathit{i}}_D(t,a)}\\ {\mathit{i}}_D(t,0)={\displaystyle {\int }_0^{\mathrm{\infty }}{\beta }_T(a)⟨\mathbf{1},{\mathit{i}}_T(t,a)⟩\mathrm{d}a{\mathbf{v}}_1}\\ {\displaystyle +{\int }_0^{\mathrm{\infty }}{\nu }_D(a)⟨{\mathbf{v}}_1,{\mathit{i}}_D(t,a)⟩\mathrm{d}a{\mathbf{v}}_2}\end{array}\end{array}$(25) ) and (26(26) $\begin{array}{r}\{\begin{array}{l}\frac{\mathrm{\partial }{\mathit{i}}_T(t,a)}{\mathrm{\partial }t}+\frac{\mathrm{\partial }{\mathit{i}}_T(t,a)}{\mathrm{\partial }a}=-\mathbf{M}{\mathit{i}}_T(t,a)\\ {\mathit{i}}_T(t,0)={\displaystyle \frac{{\overline{\mathit{S}}}_T}{{\overline{S}}_D}{\int }_0^{\mathrm{\infty }}⟨{\mathit{\beta }}_D(a),{\mathit{i}}_D(t,a)⟩\mathrm{d}a}\end{array}\end{array}$(26) ) as an abstract Cauchy problem. For this purpose, we consider the Banach spaces $X:={\mathbb{R}}^2\times L^1((0,+\mathrm{\infty }),{\mathbb{R}}^2)$, $Y={\mathbb{R}}^3\times L^1((0,+\mathrm{\infty }),{\mathbb{R}}^3)$, $X_0:=\{0_{{\mathbb{R}}^2}\}\times L^1((0,+\mathrm{\infty }),{\mathbb{R}}^2)$ and $Y_0:=\{0_{{\mathbb{R}}^3}\}\times L^1((0,+\mathrm{\infty }),{\mathbb{R}}^3)$. Let ${\mathcal{A}}_D:D({\mathcal{A}}_D)\subset X_0\to X$ and ${\mathcal{A}}_T:D({\mathcal{A}}_T)\subset Y_0\to Y$ be the linear operators defined by $\begin{array}{r}{\mathcal{A}}_D\left(\begin{array}{c}{0}_{{\mathbb{R}}^2}\\ \phi \end{array}\right)=\left(\begin{array}{c}-\phi (0)\\ -{\phi }^{\prime }-({\mu }_D+{\mathit{\theta }}_D(\cdot )\phi )\end{array}\right)\mathrm{and}{\mathcal{A}}_T\left(\begin{array}{c}{0}_{{\mathbb{R}}^3}\\ \varphi \end{array}\right)=\left(\begin{array}{c}-\varphi (0)\\ -{\varphi }^{\prime }-\mathbf{M}\varphi \end{array}\right)\end{array}$with $\begin{array}{r}D({\mathcal{A}}_D)=\{0_{{\mathbb{R}}^2}\}\times W^{1,1}((0,+\mathrm{\infty }),{\mathbb{R}}^2)\text{ and }D({\mathcal{A}}_T)=\{0_{{\mathbb{R}}^3}\}\times W^{1,1}((0,+\mathrm{\infty }),{\mathbb{R}}^3).\end{array}$Let $\mathcal{A}:D(\mathcal{A})\subset Y_0\times X_0\to Y\times X$ with $\begin{array}{r}D(\mathcal{A})=D({\mathcal{A}}_D)\times D({\mathcal{A}}_T)\mathrm{and}\mathcal{A}\left(\begin{array}{c}\mathit{\phi }\\ \mathit{\varphi }\end{array}\right)=\left(\begin{array}{c}{\mathcal{A}}_D[\mathit{\phi }]\\ {\mathcal{A}}_T[\mathit{\varphi }]\end{array}\right).\end{array}$Let $\mathcal{B}:X_0\times Y_0\to X\times Y$ be the bounded linear operator defined by (27) $\begin{array}{r}\mathcal{B}\left(\begin{array}{c}\mathit{\phi }\\ \mathit{\varphi }\end{array}\right)=\left(\begin{array}{c}\left(\begin{array}{c}{\mathbf{v}}_1{\mathcal{P}}_T[\mathit{\varphi }]+{\mathbf{v}}_2{\mathcal{Q}}_D[\mathit{\phi }]\\ 0_{L^1((0,+\mathrm{\infty }),{\mathbb{R}}^2)}\end{array}\right)\\ \frac{{\overline{\mathit{S}}}_T}{{\overline{S}}_D}\left(\begin{array}{c}{\mathcal{W}}_D[\mathit{\phi }]\\ 0_{L^1((0,+\mathrm{\infty }),{\mathbb{R}}^3)}\end{array}\right)\end{array}\right)\end{array}$(27) where we have set for each $\mathit{\phi }=(\begin{array}{c}{0}_{{\mathbb{R}}^2}\\ \phi \end{array})\in X_0$ $\begin{array}{r}\{\begin{array}{l}{\displaystyle {\mathcal{W}}_D[\mathit{\phi }]:={\int }_0^{\mathrm{\infty }}⟨{\mathit{\beta }}_D(a),\phi (a)⟩\mathrm{d}a}\\ {\displaystyle {\mathcal{P}}_D[\mathit{\phi }]:={\int }_0^{\mathrm{\infty }}\gamma (a)⟨{\mathbf{v}}_1,\phi (a)⟩\mathrm{d}a}\\ {\displaystyle {\mathcal{Q}}_D[\mathit{\phi }]:={\int }_0^{\mathrm{\infty }}{\nu }_D(a)⟨{\mathbf{v}}_1,\phi (a)⟩\mathrm{d}a}\end{array}\end{array}$and each $\mathit{\varphi }=(\begin{array}{c}{0}_{{\mathbb{R}}^3}\\ \varphi \end{array})\in Y_0$ $\begin{array}{r}\{\begin{array}{l}{\displaystyle {\mathcal{W}}_T[\mathit{\varphi }]:={\int }_0^{\mathrm{\infty }}⟨{\mathbf{e}}_3,\varphi (a)⟩\mathrm{d}a}\\ {\displaystyle {\mathcal{P}}_T[\mathit{\varphi }]:={\int }_0^{\mathrm{\infty }}{\beta }_T(a)⟨\mathbf{1},\varphi (a)⟩\mathrm{d}a.}\end{array}\end{array}$Hence, setting for all t>0 $\begin{array}{r}{\mathit{\phi }}_D(t):=\left(\begin{array}{c}{0}_{{\mathbb{R}}^2}\\ {\mathit{i}}_D(t,\cdot )\end{array}\right)\text{ and }{\mathit{\varphi }}_T(t):=\left(\begin{array}{c}{0}_{{\mathbb{R}}^3}\\ {\mathit{i}}_T(t,\cdot )\end{array}\right)\end{array}$and $\begin{array}{r}{\mathit{\phi }}_{D0}:=\left(\begin{array}{c}{0}_{{\mathbb{R}}^2}\\ {\mathit{i}}_D(0,\cdot )\end{array}\right)\text{ and }{\mathit{\varphi }}_{T0}:=\left(\begin{array}{c}{0}_{{\mathbb{R}}^3}\\ {\mathit{i}}_T(0,\cdot )\end{array}\right)\end{array}$the system (25(25) $\begin{array}{r}\{\begin{array}{l}{\displaystyle \frac{\mathrm{\partial }{\mathit{i}}_D(t,a)}{\mathrm{\partial }t}+\frac{\mathrm{\partial }{\mathit{i}}_D(t,a)}{\mathrm{\partial }a}=-({\mu }_D+\mathit{\theta }(a)){\mathit{i}}_D(t,a)}\\ {\mathit{i}}_D(t,0)={\displaystyle {\int }_0^{\mathrm{\infty }}{\beta }_T(a)⟨\mathbf{1},{\mathit{i}}_T(t,a)⟩\mathrm{d}a{\mathbf{v}}_1}\\ {\displaystyle +{\int }_0^{\mathrm{\infty }}{\nu }_D(a)⟨{\mathbf{v}}_1,{\mathit{i}}_D(t,a)⟩\mathrm{d}a{\mathbf{v}}_2}\end{array}\end{array}$(25) )–(26(26) $\begin{array}{r}\{\begin{array}{l}\frac{\mathrm{\partial }{\mathit{i}}_T(t,a)}{\mathrm{\partial }t}+\frac{\mathrm{\partial }{\mathit{i}}_T(t,a)}{\mathrm{\partial }a}=-\mathbf{M}{\mathit{i}}_T(t,a)\\ {\mathit{i}}_T(t,0)={\displaystyle \frac{{\overline{\mathit{S}}}_T}{{\overline{S}}_D}{\int }_0^{\mathrm{\infty }}⟨{\mathit{\beta }}_D(a),{\mathit{i}}_D(t,a)⟩\mathrm{d}a}\end{array}\end{array}$(26) ) can be rewritten as the following abstract Cauchy problem (28) $\begin{array}{r}\{\begin{array}{l}{\displaystyle \frac{\mathrm{d}}{\mathrm{d}t}\left(\begin{array}{c}{\mathit{\phi }}_D(t)\\ {\mathit{\varphi }}_T(t)\end{array}\right)=\mathcal{A}\left(\begin{array}{c}{\mathit{\phi }}_D(t)\\ {\mathit{\varphi }}_T(t)\end{array}\right)+\mathcal{B}\left(\begin{array}{c}{\mathit{\phi }}_D(t)\\ {\mathit{\varphi }}_T(t)\end{array}\right),t>0}\\ \left(\begin{array}{c}{\mathit{\phi }}_D(0)\\ {\mathit{\varphi }}_T(0)\end{array}\right)=\left(\begin{array}{c}{\mathit{\phi }}_{D0}\\ {\mathit{\varphi }}_{T0}\end{array}\right)\in X_0\times Y_0.\end{array}\end{array}$(28) In order to define the threshold that determines the local asymptotic stability of the disease-free equilibrium, we need first to prove the following lemma. Before proceeding, let us set $X_{+}:={\mathbb{R}}_{+}^2\times L^1((0,+\mathrm{\infty }),{\mathbb{R}}^2)$, $Y_{+}:={\mathbb{R}}_{+}^2\times L^1((0,+\mathrm{\infty }),{\mathbb{R}}^2)$, $X_{0+}:=X_{+}\cap X_0$, and $Y_{0+}:=Y_{+}\cap Y_0$.

Lemma 7.1

Let Assumption 3.1 be satisfied. Then the following properties hold true

(i)	$\mathcal{A}$ is resolvent positive i.e. $(\lambda -\mathcal{A}{)}^{-1}$ maps $X_{+}\times Y_{+}$ into itself for all large λ in the resolvent set $\rho (\mathcal{A})$ of $\mathcal{A}$;
(ii)	The spectral bound $\mathsf{s}(\mathcal{A})$ i.e. (29) $\begin{array}{r}\mathsf{s}(\mathcal{A})=sup\{\mathrm{\Re }(\lambda ):\lambda \in \sigma (\mathcal{A})\}\end{array}$(29) satisfies $\mathsf{s}(\mathcal{A})<0$ and $(\mathsf{s}(\mathcal{A}),+\mathrm{\infty })\subset \rho ({\mathcal{A}}_D)\cap \rho ({\mathcal{A}}_T)$.

Proof.

In order to obtain the conclusion of the lemma, we first prove that $\mathcal{A}$ is resolvent positive. To do this, let us note that if $\lambda \in \rho ({\mathcal{A}}_D)\cap \rho ({\mathcal{A}}_T)$ then $\lambda \in \rho (\mathcal{A})$ and $\begin{array}{r}(\lambda -\mathcal{A}{)}^{-1}\left(\begin{array}{c}\tilde{\mathit{\phi }}\\ \tilde{\mathit{\varphi }}\end{array}\right)=\left(\begin{array}{c}(\lambda -{\mathcal{A}}_D{)}^{-1}[\tilde{\mathit{\phi }}]\\ (\lambda -{\mathcal{A}}_T{)}^{-1}[\tilde{\mathit{\varphi }}]\end{array}\right)\end{array}$so that the resolvent of $\mathcal{A}$ is obtained by determining the resolvent of ${\mathcal{A}}_D$ and ${\mathcal{A}}_T$. By standard computations, we have for each $\lambda \in \rho ({\mathcal{A}}_D)$ and $\tilde{\mathit{\phi }}=(\begin{array}{c}\mathbf{x}\\ \tilde{\phi }\end{array})\in Y$ (30) $\begin{array}{r}\begin{array}{l}(\lambda -{\mathcal{A}}_D{)}^{-1}\left(\begin{array}{c}\mathbf{x}\\ \tilde{\phi }\end{array}\right)=\left(\begin{array}{c}{0}_{{\mathbb{R}}^2}\\ \phi \end{array}\right)\\ ⟺{\displaystyle \phi (a)=e^{-\lambda a}\mathbf{\Pi }(a,0)\mathbf{x}+{\int }_0^a{e}^{-\lambda (a-\tau )}\mathbf{\Pi }(a,\tau )\tilde{\phi }(\tau )\mathrm{d}\tau ,a\ge 0}\end{array}\end{array}$(30) (where $\mathbf{\Pi }(a,0)$ is defined in (9(9) $\begin{array}{r}\mathbf{\Pi }(a,\tau ):=e^{-{\int }_{\tau }^a({\mu }_D+\mathit{\theta }(l))\mathrm{d}l},\mathrm{\forall }a\ge \tau \ge 0.\end{array}$(9) )) while for each $\lambda \in \rho ({\mathcal{A}}_T)$ and $\tilde{\mathit{\varphi }}=(\begin{array}{c}\mathbf{y}\\ \tilde{\varphi }\end{array}\in X$ we have (31) $\begin{array}{r}\begin{array}{l}(\lambda -{\mathcal{A}}_T{)}^{-1}\left(\begin{array}{c}\mathbf{y}\\ \varphi \end{array}\right)=\left(\begin{array}{c}{0}_{{\mathbb{R}}^3}\\ \varphi \end{array}\right)\\ ⟺{\displaystyle \varphi (a)=e^{-\lambda a}{e}^{-\mathbf{M}a}\mathbf{y}+{\int }_0^a{e}^{-\lambda (a-\tau )}{e}^{-\mathbf{M}(a-\tau )}\tilde{\varphi }(\tau )\mathrm{d}\tau ,a\ge 0.}\end{array}\end{array}$(31) Next, recalling the definition of $\mathbf{M}$ and $\mathit{\theta }(a)$ respectively in (5(5) $\begin{array}{r}\mathbf{M}:=\left(\begin{array}{ccc}{\mu }_L+{\alpha }_L& 0& 0\\ -{\alpha }_L& {\mu }_N+{\alpha }_N& 0\\ 0& -{\alpha }_N& {\mu }_A\end{array}\right).\end{array}$(5) ) and (6(6) $\begin{array}{r}\mathit{\theta }(a):=\left(\begin{array}{cc}\gamma (a)+{\delta }_D(a)+{\nu }_D(a)& 0\\ 0& {\delta }_C(a)\end{array}\right),\end{array}$(6) ) it follows that if (32) $\begin{array}{r}\{\begin{array}{l}\lambda >-min({\mu }_L+{\alpha }_L,{\mu }_N+{\alpha }_N,{\mu }_A)=:-{\mu }_{1,-}<0\\ \mathrm{and}\lambda >-{\mu }_D-min({\mathrm{essinf}}_{{\mathbb{R}}^{+}}{\mathit{\theta }}_{11}(a),{\mathrm{essinf}}_{{\mathbb{R}}^{+}}{\mathit{\theta }}_{22}(a))=:-{\mu }_{2,-}<0\end{array}\end{array}$(32) then $\lambda \in \rho ({\mathcal{A}}_D)\cap \rho ({\mathcal{A}}_T)$. Next, recalling from [54,55] if $\mathcal{A}$ is resolvent positive then (33) $\begin{array}{r}\{\begin{array}{l}\mathsf{s}(\mathcal{A})=inf\{\lambda \in \rho (\mathcal{A}):(\lambda -\mathcal{A}{)}^{-1}(X_{+}\times Y_{+})\subseteq (X_{+}\times Y_{+})\}\\ (\mathsf{s}(\mathcal{A}),+\mathrm{\infty })\subseteq \rho (\mathcal{A})\end{array}\end{array}$(33) and we infer from (30(30) $\begin{array}{r}\begin{array}{l}(\lambda -{\mathcal{A}}_D{)}^{-1}\left(\begin{array}{c}\mathbf{x}\\ \tilde{\phi }\end{array}\right)=\left(\begin{array}{c}{0}_{{\mathbb{R}}^2}\\ \phi \end{array}\right)\\ ⟺{\displaystyle \phi (a)=e^{-\lambda a}\mathbf{\Pi }(a,0)\mathbf{x}+{\int }_0^a{e}^{-\lambda (a-\tau )}\mathbf{\Pi }(a,\tau )\tilde{\phi }(\tau )\mathrm{d}\tau ,a\ge 0}\end{array}\end{array}$(30) ) to (33(33) $\begin{array}{r}\{\begin{array}{l}\mathsf{s}(\mathcal{A})=inf\{\lambda \in \rho (\mathcal{A}):(\lambda -\mathcal{A}{)}^{-1}(X_{+}\times Y_{+})\subseteq (X_{+}\times Y_{+})\}\\ (\mathsf{s}(\mathcal{A}),+\mathrm{\infty })\subseteq \rho (\mathcal{A})\end{array}\end{array}$(33) ) that $\mathsf{s}(\mathcal{A})<0$. The proof is completed.

Thanks to Lemma 7.1 and the positiveness of the bounded linear operator $\mathcal{B}$ (i.e. it maps $X_{0+}\times Y_{0+}$ into $X_{+}\times Y_{+}$), one can use the theory developed in [56] to define a threshold which determines the sign of $\mathsf{s}(\mathcal{A}+\mathcal{B})$ i.e. the spectral bound of $\mathcal{A}+\mathcal{B}$ given by $\begin{array}{r}\mathsf{s}(\mathcal{A}+\mathcal{B})=sup\{\mathrm{\Re }(\lambda ):\lambda \in \sigma (\mathcal{A}+\mathcal{B})\}.\end{array}$More precisely, we set (34) $\begin{array}{r}{\hat{\mathcal{T}}}_0:=\mathsf{r}(\mathcal{B}(-\mathcal{A}{)}^{-1})\end{array}$(34) where $\mathsf{r}(\mathcal{B}(-\mathcal{A}{)}^{-1})$ is the spectral radius of the bounded linear operator $\mathcal{B}(-\mathcal{A}{)}^{-1}$. The following lemma will allow us to obtain the relationship between the sign of ${\hat{\mathcal{T}}}_0-1$ and the growth bound of the C${}_0$-semigroup generated by $(\mathcal{A}+\mathcal{B}{)}_0$ i.e. the part of $\mathcal{A}+\mathcal{B}$ in $X_0\times Y_0$. Note that the latter threshold ${\hat{\mathcal{T}}}_0$ does not exhibit explicitly the parameters of the model but this will be resolved after proving the local stability properties.

Lemma 7.2

Let Assumption 3.1 be satisfied. Then $\mathsf{s}(\mathcal{A}+\mathcal{B})$ and ${\hat{\mathcal{T}}}_0-1$ have the same sign. Moreover, we have the following properties

(i)	$\mathsf{s}((\mathcal{A}+\mathcal{B}{)}_0)=\mathsf{s}(\mathcal{A}+\mathcal{B})={\omega }_0((\mathcal{A}+\mathcal{B}{)}_0)$ with ${\omega }_0((\mathcal{A}+\mathcal{B}{)}_0)$ the growth bound of the C${}_0$-semigroup generated by $(\mathcal{A}+\mathcal{B}{)}_0$
(ii)	${\omega }_{0,ess}((\mathcal{A}+\mathcal{B}{)}_0)\le {\omega }_0({\mathcal{A}}_0)=\mathsf{s}({\mathcal{A}}_0)=\mathsf{s}(\mathcal{A})$ with ${\omega }_{0,ess}((\mathcal{A}+\mathcal{B}{)}_0)$ the essential growth bound of the C${}_0$-semigroup generated by $(\mathcal{A}+\mathcal{B}{)}_0$ and ${\mathcal{A}}_0$ is the part of $\mathcal{A}$ in $X_0\times Y_0$.

Proof.

We first note that $X\times Y$ is an AL-space [56, Theorem 3.14] with positive cone $X_{+}\times Y_{+}$ that is normal and generating (see [56]) so that (35) $\begin{array}{r}\mathsf{s}({\mathcal{A}}_0)={\omega }_0({\mathcal{A}}_0)\mathrm{and}\mathsf{s}((\mathcal{A}+\mathcal{B}{)}_0)={\omega }_0((\mathcal{A}+\mathcal{B}{)}_0)\end{array}$(35) with ${\omega }_0({\mathcal{A}}_0)$ the growth bound of the C${}_0$-semigroup generated by ${\mathcal{A}}_0$. The first assertion of the lemma is a direct application of [56, Theorem 3.5]. Next, we prove properties (i) and (ii). Since $\rho (\mathcal{A})=\rho ({\mathcal{A}}_0)$ and $\rho (\mathcal{A}+\mathcal{B})=\rho ((\mathcal{A}+\mathcal{B}{)}_0)$ (see [57]) it follows that $\mathsf{s}({\mathcal{A}}_0)=\mathsf{s}(\mathcal{A})$ and $\mathsf{s}((\mathcal{A}+\mathcal{B}{)}_0)=\mathsf{s}(\mathcal{A}+\mathcal{B})$. This proves (i) and the equality ${\omega }_0({\mathcal{A}}_0)=\mathsf{s}({\mathcal{A}}_0)=\mathsf{s}(\mathcal{A})$. Next, note that $\mathcal{B}$ is compact since $\mathcal{B}(X_0\times Y_0)$ is a finite-dimensional space. Therefore, we infer from [58, Theorem 1.2] that ${\omega }_{0,ess}((\mathcal{A}+\mathcal{B}{)}_0)\le {\omega }_0({\mathcal{A}}_0)$ and the proof is completed.

Lemma 7.3

Let Assumption 3.1 be satisfied. Assume in addition that the ticks' reproduction number satisfies ${\mathcal{R}}_{ST}^0>1$. Then the disease-free equilibrium is locally asymptotically stable if ${\hat{\mathcal{T}}}_0<1$ and unstable if ${\hat{\mathcal{T}}}_0>1$.

Proof.

If ${\hat{\mathcal{T}}}_0<1$ then by Lemma 7.3 we have ${\omega }_0((\mathcal{A}+\mathcal{B}{)}_0)<0$ resulting to the local asymptotic stability of the disease-free. If ${\hat{\mathcal{T}}}_0>1$ then ${\omega }_0((\mathcal{A}+\mathcal{B}{)}_0)>0$ and since ${\omega }_{0,ess}((\mathcal{A}+\mathcal{B}{)}_0)<0$ it follows that ${\omega }_{0,ess}((\mathcal{A}+\mathcal{B}{)}_0)<{\omega }_0((\mathcal{A}+\mathcal{B}{)}_0)$. Therefore, using [58, Theorem 3.2.] one knows that $(\mathcal{A}+\mathcal{B}{)}_0$ has an isolated positive eigenvalue. Thus we refer to [57, Proposition 5.7.4]) that the disease-free equilibrium is unstable.

The threshold $\widehat{{\mathcal{T}}_0}$ defined in (34(34) $\begin{array}{r}{\hat{\mathcal{T}}}_0:=\mathsf{r}(\mathcal{B}(-\mathcal{A}{)}^{-1})\end{array}$(34) ) is somehow abstract and does not allow the interpretation of the stability of the disease-free to be made in terms of the parameters of the model. Therefore, we will give in the following and equivalent threshold which incorporate explicitly the parameters of the model. To this end, we first make a remark that allows us to simplify the determination of ${\mathcal{T}}_0$. Let us set $Y_1={\mathbb{R}}^2\times L^1((0,+\mathrm{\infty }),{\mathbb{R}}^2)$, and $X_1={\mathbb{R}}^3\times L^1((0,+\mathrm{\infty }),{\mathbb{R}}^3)$. Next, we observe that $Y_1\times X_1$ is finite-dimensional and $\mathcal{B}(Y_0\times X_0)\subseteq Y_1\times X_1$ so that the spectral bound of $\mathcal{B}(-\mathcal{A}{)}^{-1}$ coincides with the spectral bound of $\mathcal{B}(-\mathcal{A}{)}_{|Y_1\times X_1}^{-1}$. The advantage of considering $\mathcal{B}(-\mathcal{A}{)}_{|Y_1\times X_1}^{-1}$ is that we can obtain the spectral radius by solving an eigenvalue problem of a matrix.

Proposition 7.4

Let Assumption 3.1 be satisfied. Assume in addition that the ticks' reproduction number satisfies ${\mathcal{R}}_{ST}^0>1$ and set (36) $\begin{array}{rl}{\mathcal{T}}_0& =\frac{1}{{\overline{S}}_D}{\int }_0^{\mathrm{\infty }}{\beta }_T(a)⟨\mathbf{1},e^{-\mathbf{M}a}{\overline{\mathit{S}}}_T⟩\mathrm{d}a[{\int }_0^{\mathrm{\infty }}{\beta }_D^1(a){\pi }_D(a)\mathrm{d}a\\ & +{\int }_0^{\mathrm{\infty }}{\nu }_D(a){\pi }_D(a)\mathrm{d}a{\int }_0^{\mathrm{\infty }}{\beta }_D^2(a){\pi }_C(a)].\end{array}$(36) Then ${\hat{\mathcal{T}}}_0-1$ and ${\mathcal{T}}_0-1$ have the same sign.

Proof.

We first give the explicit form of the linear operator $\mathcal{B}(-\mathcal{A}{)}_{|X_1\times Y_1}^{-1}$. Let $(\begin{array}{c}{\mathit{\phi }}_1\\ {\mathit{\varphi }}_1\end{array})\in X_1\times Y_1$ be given with ${\mathit{\phi }}_1:=\begin{array}{c}\mathbf{x}\\ 0_{L^1((0,+\mathrm{\infty }),{\mathbb{R}}^2)}\end{array})$ and ${\mathit{\varphi }}_1:=(\begin{array}{c}\mathbf{y}\\ 0_{L^1((0,+\mathrm{\infty }),{\mathbb{R}}^3)}\end{array})$. Then using the resolvent formula of ${\mathcal{A}}_D$ and ${\mathcal{A}}_T$ respectively in (30(30) $\begin{array}{r}\begin{array}{l}(\lambda -{\mathcal{A}}_D{)}^{-1}\left(\begin{array}{c}\mathbf{x}\\ \tilde{\phi }\end{array}\right)=\left(\begin{array}{c}{0}_{{\mathbb{R}}^2}\\ \phi \end{array}\right)\\ ⟺{\displaystyle \phi (a)=e^{-\lambda a}\mathbf{\Pi }(a,0)\mathbf{x}+{\int }_0^a{e}^{-\lambda (a-\tau )}\mathbf{\Pi }(a,\tau )\tilde{\phi }(\tau )\mathrm{d}\tau ,a\ge 0}\end{array}\end{array}$(30) ) and (31(31) $\begin{array}{r}\begin{array}{l}(\lambda -{\mathcal{A}}_T{)}^{-1}\left(\begin{array}{c}\mathbf{y}\\ \varphi \end{array}\right)=\left(\begin{array}{c}{0}_{{\mathbb{R}}^3}\\ \varphi \end{array}\right)\\ ⟺{\displaystyle \varphi (a)=e^{-\lambda a}{e}^{-\mathbf{M}a}\mathbf{y}+{\int }_0^a{e}^{-\lambda (a-\tau )}{e}^{-\mathbf{M}(a-\tau )}\tilde{\varphi }(\tau )\mathrm{d}\tau ,a\ge 0.}\end{array}\end{array}$(31) ) we obtain (37) $\begin{array}{r}(-{\mathcal{A}}_D{)}^{-1}[{\mathit{\phi }}_1]=\left(\begin{array}{c}{0}_{{\mathbb{R}}^2}\\ \mathbf{\Pi }(\cdot ,0)\mathbf{x}\end{array}\right)=:{\mathit{\phi }}_0\end{array}$(37) and (38) $\begin{array}{r}(-{\mathcal{A}}_T{)}^{-1}[{\mathit{\varphi }}_1]=\left(\begin{array}{c}{0}_{{\mathbb{R}}^3}\\ e^{-\mathbf{M}\cdot }\mathbf{y}\end{array}\right)=:{\mathit{\varphi }}_0.\end{array}$(38) Hence using the explicit form of $\mathcal{B}$ in (27(27) $\begin{array}{r}\mathcal{B}\left(\begin{array}{c}\mathit{\phi }\\ \mathit{\varphi }\end{array}\right)=\left(\begin{array}{c}\left(\begin{array}{c}{\mathbf{v}}_1{\mathcal{P}}_T[\mathit{\varphi }]+{\mathbf{v}}_2{\mathcal{Q}}_D[\mathit{\phi }]\\ 0_{L^1((0,+\mathrm{\infty }),{\mathbb{R}}^2)}\end{array}\right)\\ \frac{{\overline{\mathit{S}}}_T}{{\overline{S}}_D}\left(\begin{array}{c}{\mathcal{W}}_D[\mathit{\phi }]\\ 0_{L^1((0,+\mathrm{\infty }),{\mathbb{R}}^3)}\end{array}\right)\end{array}\right)\end{array}$(27) ) it comes (39) $\begin{array}{r}\mathcal{B}(-\mathcal{A}{)}^{-1}\left(\begin{array}{c}{\mathit{\phi }}_1\\ {\mathit{\varphi }}_1\end{array}\right)=\mathcal{B}\left(\begin{array}{c}{\mathit{\phi }}_0\\ {\mathit{\varphi }}_0\end{array}\right)=\left(\begin{array}{c}\left(\begin{array}{c}{\mathbf{v}}_1{\mathcal{P}}_T[{\mathit{\varphi }}_0]+{\mathbf{v}}_2{\mathcal{Q}}_D[{\mathit{\phi }}_0]\\ 0_{L^1((0,+\mathrm{\infty }),{\mathbb{R}}^2)}\end{array}\right)\\ {\displaystyle \frac{{\overline{\mathit{S}}}_T}{{\overline{S}}_D}\left(\begin{array}{c}{\mathcal{W}}_D[{\mathit{\phi }}_0]\\ 0_{L^1((0,+\mathrm{\infty }),{\mathbb{R}}^3)}\end{array}\right)}\end{array}\right)\end{array}$(39) with $\begin{array}{r}\{\begin{array}{l}{\mathbf{v}}_1{\mathcal{P}}_T[{\mathit{\varphi }}_0]+{\mathbf{v}}_2{\mathcal{Q}}_D[{\mathit{\phi }}_0]={\mathbf{v}}_1{\int }_0^{\mathrm{\infty }}{\beta }_T(a)⟨\mathbf{1},e^{-\mathbf{M}a}\mathbf{y}⟩\mathrm{d}a+{\mathbf{v}}_2{\int }_0^{\mathrm{\infty }}{\nu }_D(a)⟨{\mathbf{v}}_1,\mathbf{\Pi }(a,0)\mathbf{x}⟩\mathrm{d}a\\ {\mathcal{W}}_D[{\mathit{\phi }}_0]={\int }_0^{\mathrm{\infty }}⟨{\mathit{\beta }}_D(a),\mathbf{\Pi }(a,0)\mathbf{x}⟩\mathrm{d}a\end{array}\end{array}$so that the spectral radius of $\mathcal{B}(-\mathcal{A}{)}^{-1}$ is given by the spectral radius of the linear operator $\mathcal{C}:{\mathbb{R}}^2\times {\mathbb{R}}^3\to {\mathbb{R}}^2\times {\mathbb{R}}^3$ given by $\begin{array}{r}\mathcal{C}\left(\begin{array}{c}\mathbf{x}\\ \mathbf{y}\end{array}\right)=\left(\begin{array}{c}{\displaystyle {\mathbf{v}}_1{\int }_0^{\mathrm{\infty }}{\beta }_T(a)⟨\mathbf{1},e^{-\mathbf{M}a}\mathbf{y}⟩\mathrm{d}a+{\mathbf{v}}_2{\int }_0^{\mathrm{\infty }}{\nu }_D(a)⟨{\mathbf{v}}_1,\mathbf{\Pi }(a,0)\mathbf{x}⟩\mathrm{d}a}\\ {\displaystyle \frac{{\overline{\mathit{S}}}_T}{{\overline{S}}_D}{\int }_0^{\mathrm{\infty }}⟨{\mathit{\beta }}_D(a),\mathbf{\Pi }(a,0)\mathbf{x}⟩\mathrm{d}a}\end{array}\right)\end{array}$Since $\mathcal{C}$ is a positive linear operator of finite-dimensional spaces, its spectral radius ${\hat{\mathcal{T}}}_0$ is an eigenvalue. Moreover, let $(\begin{array}{c}\mathbf{x}\\ \mathbf{y}\end{array})\in {\mathbb{R}}^2\times {\mathbb{R}}^3$ be a non zero vector such that $\mathcal{C}(\begin{array}{c}\mathbf{x}\\ \mathbf{y}\end{array})={\hat{\mathcal{T}}}_0(\begin{array}{c}\mathbf{x}\\ \mathbf{y})\end{array}$. Then we have the following system of equations (40) $\begin{array}{r}\{\begin{array}{l}{\displaystyle {\mathbf{v}}_1{\int }_0^{\mathrm{\infty }}{\beta }_T(a)⟨\mathbf{1},e^{-\mathbf{M}a}\mathbf{y}⟩\mathrm{d}a+{\mathbf{v}}_2{\int }_0^{\mathrm{\infty }}{\nu }_D(a)⟨{\mathbf{v}}_1,\mathbf{\Pi }(a,0)\mathbf{x}⟩\mathrm{d}a={\hat{\mathcal{T}}}_0\mathbf{x}}\\ {\displaystyle \frac{{\overline{\mathit{S}}}_T}{{\overline{S}}_D}{\int }_0^{\mathrm{\infty }}⟨{\mathit{\beta }}_D(a),\mathbf{\Pi }(a,0)\mathbf{x}⟩\mathrm{d}a={\hat{\mathcal{T}}}_0\mathbf{y}}\end{array}\end{array}$(40) from where we obtain the equality (41) $\begin{array}{r}\mathbf{y}=\frac{1}{{\hat{\mathcal{T}}}_0}\frac{{\overline{\mathit{S}}}_T}{{\overline{S}}_D}{\int }_0^{\mathrm{\infty }}⟨{\mathit{\beta }}_D(a),\mathbf{\Pi }(a,0)\mathbf{x}⟩\mathrm{d}a.\end{array}$(41) Hence, plugging (41(41) $\begin{array}{r}\mathbf{y}=\frac{1}{{\hat{\mathcal{T}}}_0}\frac{{\overline{\mathit{S}}}_T}{{\overline{S}}_D}{\int }_0^{\mathrm{\infty }}⟨{\mathit{\beta }}_D(a),\mathbf{\Pi }(a,0)\mathbf{x}⟩\mathrm{d}a.\end{array}$(41) ) into the first equation of (40(40) $\begin{array}{r}\{\begin{array}{l}{\displaystyle {\mathbf{v}}_1{\int }_0^{\mathrm{\infty }}{\beta }_T(a)⟨\mathbf{1},e^{-\mathbf{M}a}\mathbf{y}⟩\mathrm{d}a+{\mathbf{v}}_2{\int }_0^{\mathrm{\infty }}{\nu }_D(a)⟨{\mathbf{v}}_1,\mathbf{\Pi }(a,0)\mathbf{x}⟩\mathrm{d}a={\hat{\mathcal{T}}}_0\mathbf{x}}\\ {\displaystyle \frac{{\overline{\mathit{S}}}_T}{{\overline{S}}_D}{\int }_0^{\mathrm{\infty }}⟨{\mathit{\beta }}_D(a),\mathbf{\Pi }(a,0)\mathbf{x}⟩\mathrm{d}a={\hat{\mathcal{T}}}_0\mathbf{y}}\end{array}\end{array}$(40) ) gives $\begin{array}{rl}& {\mathbf{v}}_1\frac{1}{{\hat{\mathcal{T}}}_0}\frac{1}{{\overline{S}}_D}{\int }_0^{\mathrm{\infty }}{\beta }_T(a)⟨\mathbf{1},e^{-\mathbf{M}a}{\overline{\mathit{S}}}_T⟩\mathrm{d}a{\int }_0^{\mathrm{\infty }}⟨{\mathit{\beta }}_D(a),\mathbf{\Pi }(a,0)\mathbf{x}⟩\mathrm{d}a\\ & +{\mathbf{v}}_2{\int }_0^{\mathrm{\infty }}{\nu }_D(a)⟨{\mathbf{v}}_1,\mathbf{\Pi }(a,0)\mathbf{x}⟩\mathrm{d}a={\hat{\mathcal{T}}}_0\mathbf{x}.\end{array}$Thus, setting $\mathbf{x}=x_1{\mathbf{v}}_1+x_2{\mathbf{v}}_2$ and recalling that $⟨{\mathit{\beta }}_D(a),\mathbf{\Pi }(a,0){\mathbf{v}}_1⟩={\beta }_D^1(a){\pi }_D(a)$, $⟨{\mathit{\beta }}_D(a),\mathbf{\Pi }(a,0){\mathbf{v}}_2⟩={\beta }_D^2(a){\pi }_C(a)$, and $⟨{\mathbf{v}}_1,\mathbf{\Pi }(a,0){\mathbf{v}}_1⟩={\pi }_D(a)$ we obtain by identification $\begin{array}{r}\{\begin{array}{l}{\displaystyle \frac{1}{{\overline{S}}_D}{\int }_0^{\mathrm{\infty }}{\beta }_T(a)⟨\mathbf{1},e^{-\mathbf{M}a}{\overline{\mathit{S}}}_T⟩\mathrm{d}a{\int }_0^{\mathrm{\infty }}[x_1{\beta }_D^1(a){\pi }_D(a)+x_2{\beta }_D^2(a){\pi }_C(a)]\mathrm{d}a={\hat{\mathcal{T}}}_0^2{x}_1}\\ {\displaystyle x_1{\int }_0^{\mathrm{\infty }}{\nu }_D(a){\pi }_D(a)\mathrm{d}a={\hat{\mathcal{T}}}_0{x}_2}\end{array}\end{array}$so that (42) $\begin{array}{rl}& \frac{1}{{\overline{S}}_D}{\int }_0^{\mathrm{\infty }}{\beta }_T(a)⟨\mathbf{1},e^{-\mathbf{M}a}{\overline{\mathit{S}}}_T⟩\mathrm{d}a[x_1{\int }_0^{\mathrm{\infty }}{\beta }_D^1(a){\pi }_D(a)\mathrm{d}a\\ & +\frac{1}{{\hat{\mathcal{T}}}_0}{x}_1{\int }_0^{\mathrm{\infty }}{\nu }_D(a){\pi }_D(a)\mathrm{d}a{\int }_0^{\mathrm{\infty }}{\beta }_D^2(a){\pi }_C(a)\mathrm{d}a]={\hat{\mathcal{T}}}_0^2{x}_1.\end{array}$(42) From the above computations, we see that $(\begin{array}{c}\mathbf{x}\\ \mathbf{y}\end{array})$ is the null vector if and only if $x_1=0$. Therefore, the spectral radius of $\mathcal{C}$ satisfies (43) $\begin{array}{r}{\hat{\mathcal{T}}}_0^2=\frac{1}{{\hat{\mathcal{T}}}_0}{\mathcal{T}}_{0,2}+{\mathcal{T}}_{0,3}\end{array}$(43) with $\begin{array}{r}\{\begin{array}{l}{\displaystyle {\mathcal{T}}_{0,2}:=\frac{1}{{\overline{S}}_D}{\int }_0^{\mathrm{\infty }}{\beta }_T(a)⟨\mathbf{1},e^{-\mathbf{M}a}{\overline{\mathit{S}}}_T⟩\mathrm{d}a{\int }_0^{\mathrm{\infty }}{\beta }_D^1(a){\pi }_D(a)\mathrm{d}a}\\ {\displaystyle {\mathcal{T}}_{0,3}:=\frac{1}{{\overline{S}}_D}{\int }_0^{\mathrm{\infty }}{\beta }_T(a)⟨\mathbf{1},e^{-\mathbf{M}a}{\overline{\mathit{S}}}_T⟩\mathrm{d}a{\int }_0^{\mathrm{\infty }}{\nu }_D(a){\pi }_D(a)\mathrm{d}a{\int }_0^{\mathrm{\infty }}{\beta }_D^2(a){\pi }_C(a)\mathrm{d}a.}\end{array}\end{array}$Next, observe that ${\mathcal{T}}_0$ defined in (15(15) $\begin{array}{r}{\mathcal{T}}_0={\mathcal{T}}_0^{D\to T}\times {\mathcal{T}}_0^{T\to D},\end{array}$(15) ) is also given by ${\mathcal{T}}_0={\mathcal{T}}_{0,2}+{\mathcal{T}}_{0,3}$. Moreover, by similar computations, one obtains that the non-zero eigenvalues of $\mathcal{C}$ satisfy (44) $\begin{array}{r}{\lambda }^3={\mathcal{T}}_{0,2}+\lambda {\mathcal{T}}_{0,3}.\end{array}$(44) To prove that ${\hat{\mathcal{T}}}_0-1$ and ${\mathcal{T}}_0-1$ have the same sign, we will make use of (43(43) $\begin{array}{r}{\hat{\mathcal{T}}}_0^2=\frac{1}{{\hat{\mathcal{T}}}_0}{\mathcal{T}}_{0,2}+{\mathcal{T}}_{0,3}\end{array}$(43) ) and (44(44) $\begin{array}{r}{\lambda }^3={\mathcal{T}}_{0,2}+\lambda {\mathcal{T}}_{0,3}.\end{array}$(44) ). Indeed if ${\hat{\mathcal{T}}}_0=1$ then (43(43) $\begin{array}{r}{\hat{\mathcal{T}}}_0^2=\frac{1}{{\hat{\mathcal{T}}}_0}{\mathcal{T}}_{0,2}+{\mathcal{T}}_{0,3}\end{array}$(43) ) implies that ${\mathcal{T}}_0=1$. If ${\mathcal{T}}_0=1$ then (44(44) $\begin{array}{r}{\lambda }^3={\mathcal{T}}_{0,2}+\lambda {\mathcal{T}}_{0,3}.\end{array}$(44) ) becomes ${\lambda }^3=(1-{\mathcal{T}}_{0,3})+\lambda {\mathcal{T}}_{0,3}$ with $1-{\mathcal{T}}_{0,3}\ge 0$ so that the non-zero solution of (44(44) $\begin{array}{r}{\lambda }^3={\mathcal{T}}_{0,2}+\lambda {\mathcal{T}}_{0,3}.\end{array}$(44) ) have modulus less than one. This proves that ${\hat{\mathcal{T}}}_0=1$ whenever ${\mathcal{T}}_0=1$. To complete the proof, we observe that if ${\mathcal{T}}_0>1$ then (43(43) $\begin{array}{r}{\hat{\mathcal{T}}}_0^2=\frac{1}{{\hat{\mathcal{T}}}_0}{\mathcal{T}}_{0,2}+{\mathcal{T}}_{0,3}\end{array}$(43) ) implies that $1<{\mathcal{T}}_0^2=\frac{1}{{\mathcal{T}}_0}{\mathcal{T}}_{0,2}+{\mathcal{T}}_{0,3}<{\mathcal{T}}_{0,2}+{\mathcal{T}}_{0,3}={\mathcal{R}}_0$. Similarly if ${\mathcal{T}}_0<1$ then ${\mathcal{R}}_0<1$.

8. Proof of Theorem 3.4

In order to determine the endemic equilibrium of the model (7(7) $\begin{array}{r}\{\begin{array}{l}{\displaystyle \frac{\mathrm{d}{S}_D(t)}{\mathrm{d}t}={\theta }_D-\frac{S_D(t)}{N_D(t)}{\int }_0^{\mathrm{\infty }}{\beta }_T(a)⟨\mathbf{1},{\mathit{i}}_T(t,a)⟩\mathrm{d}a-{\mu }_D{S}_D(t),}\\ {\displaystyle \frac{\mathrm{\partial }{\mathit{i}}_D(t,a)}{\mathrm{\partial }t}+\frac{\mathrm{\partial }{\mathit{i}}_D(t,a)}{\mathrm{\partial }a}=-({\mu }_D+\mathit{\theta }(a)){\mathit{i}}_D(t,a),}\\ {\mathit{i}}_D(t,0)={\displaystyle \frac{S_D(t)+ϵR_D(t)}{N_D(t)}{\int }_0^{\mathrm{\infty }}{\beta }_T(a)⟨\mathbf{1},{\mathit{i}}_T(t,a)⟩\mathrm{d}a{\mathbf{v}}_1}\\ {\displaystyle +{\int }_0^{\mathrm{\infty }}{\nu }_D(a)⟨{\mathbf{v}}_1,{\mathit{i}}_D(t,a)⟩\mathrm{d}a{\mathbf{v}}_2,}\\ {\displaystyle \frac{\mathrm{d}{R}_D(t)}{\mathrm{d}t}={\int }_0^{\mathrm{\infty }}\gamma (a)⟨{\mathbf{v}}_1,{\mathit{i}}_D(t,a)⟩\mathrm{d}a-ϵ\frac{R_D(t)}{N_D(t)}{\int }_0^{\mathrm{\infty }}{\beta }_T(a)⟨\mathbf{1},{\mathit{i}}_T(t,a)⟩\mathrm{d}a-{\mu }_D{R}_D(t),}\end{array}\end{array}$(7) )–(8(8) $\begin{array}{r}\{\begin{array}{l}{\displaystyle \frac{\mathrm{d}{S}_{TE}(t)}{\mathrm{d}t}=r_E(1-\frac{S_{TE}(t)}{K}){\displaystyle (⟨{\mathit{S}}_T(t),{\mathbf{e}}_3⟩+{\int }_0^{\mathrm{\infty }}⟨{\mathit{i}}_T(t,a),{\mathbf{e}}_3⟩\mathrm{d}a)-({\alpha }_E+{\mu }_E)S_{TE}(t),}}\\ {\displaystyle \frac{\mathrm{d}{\mathit{S}}_T(t)}{\mathrm{d}t}={\alpha }_E{S}_{TE}(t){\mathbf{e}}_1-\mathbf{M}{\mathit{S}}_T(t)-\frac{{\mathit{S}}_T(t)}{N_D(t)}{\int }_0^{\mathrm{\infty }}⟨{\mathit{\beta }}_D(a),{\mathit{i}}_D(t,a)⟩\mathrm{d}a,}\\ {\displaystyle \frac{\mathrm{\partial }{\mathit{i}}_T(t,a)}{\mathrm{\partial }t}+\frac{\mathrm{\partial }{\mathit{i}}_T(t,a)}{\mathrm{\partial }a}=-\mathbf{M}{\mathit{i}}_T(t,a),}\\ {\mathit{i}}_T(t,0)={\displaystyle \frac{{\mathit{S}}_T(t)}{N_D(t)}{\int }_0^{\mathrm{\infty }}⟨{\mathit{\beta }}_D(a),{\mathit{i}}_D(t,a)⟩\mathrm{d}a,}\end{array}\end{array}$(8) ), we proceed in several steps. In the first step, we give an equivalent system to (45(45) $\begin{array}{r}\{\begin{array}{l}{\displaystyle 0=r_E(1-\frac{S_{TE}^{\ast }}{K}){\displaystyle (⟨{\mathit{S}}_T^{\ast },{\mathbf{e}}_3⟩+{\int }_0^{\mathrm{\infty }}⟨{\mathit{i}}_T^{\ast }(a),{\mathbf{e}}_3⟩\mathrm{d}a)-({\alpha }_E+{\mu }_E)S_{TE}^{\ast }}}\\ {\displaystyle {\mathit{i}}_T^{\ast }(a)=\frac{1}{N_D^{\ast }}{\int }_0^{\mathrm{\infty }}⟨{\mathit{\beta }}_D(a),{\mathit{i}}_D^{\ast }(a)⟩\mathrm{d}a{e}^{-\mathbf{M}a}{\mathit{S}}_T^{\ast }}\\ {\displaystyle 0_{{\mathbb{R}}^3}={\alpha }_E{S}_{TE}^{\ast }{\mathbf{e}}_1-\frac{1}{N_D^{\ast }}{\int }_0^{\mathrm{\infty }}⟨{\mathit{\beta }}_D(a),{\mathit{i}}_D^{\ast }(a)⟩\mathrm{d}a{\mathit{S}}_T^{\ast }-\mathbf{M}{\mathit{S}}_T^{,}}\end{array}\end{array}$(45) ), the second step is concerned with the solvability of (46(46) $\begin{array}{r}\{\begin{array}{l}{\displaystyle 0={\theta }_D-\frac{S_D^{\ast }}{N_D^{\ast }}{\int }_0^{\mathrm{\infty }}{\beta }_T(a)⟨\mathbf{1},{\mathit{i}}_T^{\ast }(a)⟩\mathrm{d}a-{\mu }_D{S}_D^{\ast }}\\ {\displaystyle {\mathit{i}}_D^{\ast }(a)=\frac{S_D^{\ast }+ϵR_D^{\ast }}{N_D^{\ast }}{\int }_0^{\mathrm{\infty }}{\beta }_T(a)⟨\mathbf{1},{\mathit{i}}_T^{\ast }(a)⟩\mathrm{d}a{\pi }_D(a){\mathbf{v}}_1}\\ {\displaystyle +{\int }_0^{\mathrm{\infty }}{\nu }_D(\tau )⟨{\mathbf{v}}_1,{\mathit{i}}_D^{\ast }(a)⟩\mathrm{d}a{\pi }_C(a){\mathbf{v}}_2}\\ {\displaystyle 0={\int }_0^{\mathrm{\infty }}\gamma (a)⟨{\mathbf{v}}_1,{\mathit{i}}_D^{\ast }(a)⟩\mathrm{d}a-ϵ\frac{R_D^{\ast }}{N_D^{\ast }}{\int }_0^{\mathrm{\infty }}{\beta }_T(a)⟨\mathbf{1},{\mathit{i}}_T^{\ast }(a)⟩\mathrm{d}a-{\mu }_D{R}_D^{\ast }}\\ {\displaystyle 0={\theta }_D-{\mu }_D{N}_D^{\ast }-{\int }_0^{\mathrm{\infty }}{\delta }_D(a)⟨{\mathbf{v}}_1,{\mathit{i}}_D^{\ast }(a)⟩\mathrm{d}a-{\int }_0^{\mathrm{\infty }}{\delta }_C(a)⟨{\mathbf{v}}_2,{\mathit{i}}_D^{\ast }(a)⟩\mathrm{d}a.}\end{array}\end{array}$(46) ), and the last step concerns the bifurcation properties. In this section, we will always suppose that Assumption 3.1 is satisfied and ${\mathcal{R}}_{ST}^0>1$. Let us note that using the notation (18(18) $\begin{array}{r}\{\begin{array}{l}{\pi }_D(a)=e^{-{\int }_0^a({\mu }_D+{\delta }_D(l)+{\nu }_D(l))\mathrm{d}l}=\mathbf{\Pi }(a,0){\mathbf{v}}_1,\\ {\pi }_C(a)=e^{-{\int }_0^a({\mu }_D+{\delta }_C(l))\mathrm{d}l}=\mathbf{\Pi }(a,0){\mathbf{v}}_2.\end{array}\end{array}$(18) ) together with Theorem 3.2, one knows that an endemic equilibrium ${\mathcal{E}}^{\ast }$ defined in (19(19) $\begin{array}{r}{\mathcal{E}}^{\ast }=(S_{TE}^{\ast },{\mathit{S}}_T^{\ast },{\mathit{i}}_T^{\ast },S_D^{\ast },{\mathit{i}}_D^{\ast },R_D^{\ast })\end{array}$(19) ) satisfies (45) $\begin{array}{r}\{\begin{array}{l}{\displaystyle 0=r_E(1-\frac{S_{TE}^{\ast }}{K}){\displaystyle (⟨{\mathit{S}}_T^{\ast },{\mathbf{e}}_3⟩+{\int }_0^{\mathrm{\infty }}⟨{\mathit{i}}_T^{\ast }(a),{\mathbf{e}}_3⟩\mathrm{d}a)-({\alpha }_E+{\mu }_E)S_{TE}^{\ast }}}\\ {\displaystyle {\mathit{i}}_T^{\ast }(a)=\frac{1}{N_D^{\ast }}{\int }_0^{\mathrm{\infty }}⟨{\mathit{\beta }}_D(a),{\mathit{i}}_D^{\ast }(a)⟩\mathrm{d}a{e}^{-\mathbf{M}a}{\mathit{S}}_T^{\ast }}\\ {\displaystyle 0_{{\mathbb{R}}^3}={\alpha }_E{S}_{TE}^{\ast }{\mathbf{e}}_1-\frac{1}{N_D^{\ast }}{\int }_0^{\mathrm{\infty }}⟨{\mathit{\beta }}_D(a),{\mathit{i}}_D^{\ast }(a)⟩\mathrm{d}a{\mathit{S}}_T^{\ast }-\mathbf{M}{\mathit{S}}_T^{,}}\end{array}\end{array}$(45) and (46) $\begin{array}{r}\{\begin{array}{l}{\displaystyle 0={\theta }_D-\frac{S_D^{\ast }}{N_D^{\ast }}{\int }_0^{\mathrm{\infty }}{\beta }_T(a)⟨\mathbf{1},{\mathit{i}}_T^{\ast }(a)⟩\mathrm{d}a-{\mu }_D{S}_D^{\ast }}\\ {\displaystyle {\mathit{i}}_D^{\ast }(a)=\frac{S_D^{\ast }+ϵR_D^{\ast }}{N_D^{\ast }}{\int }_0^{\mathrm{\infty }}{\beta }_T(a)⟨\mathbf{1},{\mathit{i}}_T^{\ast }(a)⟩\mathrm{d}a{\pi }_D(a){\mathbf{v}}_1}\\ {\displaystyle +{\int }_0^{\mathrm{\infty }}{\nu }_D(\tau )⟨{\mathbf{v}}_1,{\mathit{i}}_D^{\ast }(a)⟩\mathrm{d}a{\pi }_C(a){\mathbf{v}}_2}\\ {\displaystyle 0={\int }_0^{\mathrm{\infty }}\gamma (a)⟨{\mathbf{v}}_1,{\mathit{i}}_D^{\ast }(a)⟩\mathrm{d}a-ϵ\frac{R_D^{\ast }}{N_D^{\ast }}{\int }_0^{\mathrm{\infty }}{\beta }_T(a)⟨\mathbf{1},{\mathit{i}}_T^{\ast }(a)⟩\mathrm{d}a-{\mu }_D{R}_D^{\ast }}\\ {\displaystyle 0={\theta }_D-{\mu }_D{N}_D^{\ast }-{\int }_0^{\mathrm{\infty }}{\delta }_D(a)⟨{\mathbf{v}}_1,{\mathit{i}}_D^{\ast }(a)⟩\mathrm{d}a-{\int }_0^{\mathrm{\infty }}{\delta }_C(a)⟨{\mathbf{v}}_2,{\mathit{i}}_D^{\ast }(a)⟩\mathrm{d}a.}\end{array}\end{array}$(46) On the stationary states of ticks: Let us first prove the following lemma which describes the stationary states of the eggs.

Lemma 8.1

The total number of eggs remains unchanged at the endemic and the disease-free equilibrium i.e. ${\overline{S}}_{TE}^{\ast }={\overline{S}}_{TE}$. Moreover, we always have (47) $\begin{array}{r}{\mathit{S}}_T^{\ast }+{\int }_0^{\mathrm{\infty }}{\mathit{i}}_T^{\ast }(a)\mathrm{d}a={\overline{\mathit{S}}}_T.\end{array}$(47)

Proof.

Note that, using the third equation of (45(45) $\begin{array}{r}\{\begin{array}{l}{\displaystyle 0=r_E(1-\frac{S_{TE}^{\ast }}{K}){\displaystyle (⟨{\mathit{S}}_T^{\ast },{\mathbf{e}}_3⟩+{\int }_0^{\mathrm{\infty }}⟨{\mathit{i}}_T^{\ast }(a),{\mathbf{e}}_3⟩\mathrm{d}a)-({\alpha }_E+{\mu }_E)S_{TE}^{\ast }}}\\ {\displaystyle {\mathit{i}}_T^{\ast }(a)=\frac{1}{N_D^{\ast }}{\int }_0^{\mathrm{\infty }}⟨{\mathit{\beta }}_D(a),{\mathit{i}}_D^{\ast }(a)⟩\mathrm{d}a{e}^{-\mathbf{M}a}{\mathit{S}}_T^{\ast }}\\ {\displaystyle 0_{{\mathbb{R}}^3}={\alpha }_E{S}_{TE}^{\ast }{\mathbf{e}}_1-\frac{1}{N_D^{\ast }}{\int }_0^{\mathrm{\infty }}⟨{\mathit{\beta }}_D(a),{\mathit{i}}_D^{\ast }(a)⟩\mathrm{d}a{\mathit{S}}_T^{\ast }-\mathbf{M}{\mathit{S}}_T^{,}}\end{array}\end{array}$(45) ) we have (48) $\begin{array}{r}{0}_{{\mathbb{R}}^3}={\alpha }_E{S}_{TE}{\mathbf{M}}^{-1}{\mathbf{e}}_1-\frac{1}{N_D^{\ast }}{\int }_0^{\mathrm{\infty }}⟨{\mathit{\beta }}_D(a),{\mathit{i}}_D^{\ast }(a)⟩\mathrm{d}a{\mathbf{M}}^{-1}{\mathit{S}}_T^{\ast }-{\mathit{S}}_T^{\ast }\end{array}$(48) Next, integrating the second equation of (45(45) $\begin{array}{r}\{\begin{array}{l}{\displaystyle 0=r_E(1-\frac{S_{TE}^{\ast }}{K}){\displaystyle (⟨{\mathit{S}}_T^{\ast },{\mathbf{e}}_3⟩+{\int }_0^{\mathrm{\infty }}⟨{\mathit{i}}_T^{\ast }(a),{\mathbf{e}}_3⟩\mathrm{d}a)-({\alpha }_E+{\mu }_E)S_{TE}^{\ast }}}\\ {\displaystyle {\mathit{i}}_T^{\ast }(a)=\frac{1}{N_D^{\ast }}{\int }_0^{\mathrm{\infty }}⟨{\mathit{\beta }}_D(a),{\mathit{i}}_D^{\ast }(a)⟩\mathrm{d}a{e}^{-\mathbf{M}a}{\mathit{S}}_T^{\ast }}\\ {\displaystyle 0_{{\mathbb{R}}^3}={\alpha }_E{S}_{TE}^{\ast }{\mathbf{e}}_1-\frac{1}{N_D^{\ast }}{\int }_0^{\mathrm{\infty }}⟨{\mathit{\beta }}_D(a),{\mathit{i}}_D^{\ast }(a)⟩\mathrm{d}a{\mathit{S}}_T^{\ast }-\mathbf{M}{\mathit{S}}_T^{,}}\end{array}\end{array}$(45) ) it follows that (49) $\begin{array}{rl}{\int }_0^{\mathrm{\infty }}{\mathit{i}}_T^{\ast }(a)\mathrm{d}a& =\frac{1}{N_D^{\ast }}{\int }_0^{\mathrm{\infty }}⟨{\mathit{\beta }}_D(a),{\mathit{i}}_D^{\ast }(a)⟩\mathrm{d}a{\int }_0^{\mathrm{\infty }}{e}^{-\mathbf{M}a}{\mathit{S}}_T^{\ast }\mathrm{d}a\\ & =\frac{1}{N_D^{\ast }}{\int }_0^{\mathrm{\infty }}⟨{\mathit{\beta }}_D(a),{\mathit{i}}_D^{\ast }(a)⟩\mathrm{d}a{\mathbf{M}}^{-1}{\mathit{S}}_T^{\ast }\end{array}$(49) Thus, by combining (48(48) $\begin{array}{r}{0}_{{\mathbb{R}}^3}={\alpha }_E{S}_{TE}{\mathbf{M}}^{-1}{\mathbf{e}}_1-\frac{1}{N_D^{\ast }}{\int }_0^{\mathrm{\infty }}⟨{\mathit{\beta }}_D(a),{\mathit{i}}_D^{\ast }(a)⟩\mathrm{d}a{\mathbf{M}}^{-1}{\mathit{S}}_T^{\ast }-{\mathit{S}}_T^{\ast }\end{array}$(48) ) and (49(49) $\begin{array}{rl}{\int }_0^{\mathrm{\infty }}{\mathit{i}}_T^{\ast }(a)\mathrm{d}a& =\frac{1}{N_D^{\ast }}{\int }_0^{\mathrm{\infty }}⟨{\mathit{\beta }}_D(a),{\mathit{i}}_D^{\ast }(a)⟩\mathrm{d}a{\int }_0^{\mathrm{\infty }}{e}^{-\mathbf{M}a}{\mathit{S}}_T^{\ast }\mathrm{d}a\\ & =\frac{1}{N_D^{\ast }}{\int }_0^{\mathrm{\infty }}⟨{\mathit{\beta }}_D(a),{\mathit{i}}_D^{\ast }(a)⟩\mathrm{d}a{\mathbf{M}}^{-1}{\mathit{S}}_T^{\ast }\end{array}$(49) ) we obtain (50) $\begin{array}{r}{\int }_0^{\mathrm{\infty }}{\mathit{i}}_T^{\ast }(a)\mathrm{d}a+{\mathit{S}}_T^{\ast }={\alpha }_E{S}_{TE}^{\ast }{\mathbf{M}}^{-1}{\mathbf{e}}_1=\frac{S_{TE}^{\ast }}{{\overline{S}}_{TE}}{\overline{\mathit{S}}}_T.\end{array}$(50) Therefore, plugging (50(50) $\begin{array}{r}{\int }_0^{\mathrm{\infty }}{\mathit{i}}_T^{\ast }(a)\mathrm{d}a+{\mathit{S}}_T^{\ast }={\alpha }_E{S}_{TE}^{\ast }{\mathbf{M}}^{-1}{\mathbf{e}}_1=\frac{S_{TE}^{\ast }}{{\overline{S}}_{TE}}{\overline{\mathit{S}}}_T.\end{array}$(50) ) in the first equation of (53(53) $\begin{array}{r}\{\begin{array}{l}{\displaystyle {\mathit{i}}_T^{\ast }(a)=\frac{1}{N_D^{\ast }}{\int }_0^{\mathrm{\infty }}⟨{\mathit{\beta }}_D(a),{\mathit{i}}_D^{\ast }(a)⟩\mathrm{d}a{e}^{-\mathbf{M}a}{\mathit{S}}_T^{\ast }}\\ {\displaystyle 0_{{\mathbb{R}}^3}={\alpha }_E{S}_{TE}{\mathbf{e}}_1-\frac{1}{N_D^{\ast }}{\int }_0^{\mathrm{\infty }}⟨{\mathit{\beta }}_D(a),{\mathit{i}}_D^{\ast }(a)⟩\mathrm{d}a{\mathit{S}}_T^{\ast }-\mathbf{M}{\mathit{S}}_T^{\ast }.}\end{array}\end{array}$(53) ) we obtain the following equality (51) $\begin{array}{r}0=r_E(1-\frac{S_{TE}^{\ast }}{K}){\displaystyle \frac{S_{TE}^{\ast }}{{\overline{S}}_{TE}}⟨{\overline{\mathit{S}}}_T,{\mathbf{e}}_3⟩-({\alpha }_E+{\mu }_E)S_{TE}^{\ast }}\end{array}$(51) or equivalently (52) $\begin{array}{r}0=r_E(1-\frac{S_{TE}^{\ast }}{K}){\displaystyle ⟨{\overline{\mathit{S}}}_T,{\mathbf{e}}_3⟩-({\alpha }_E+{\mu }_E){\overline{S}}_{TE}.}\end{array}$(52) The right-hand side of (52(52) $\begin{array}{r}0=r_E(1-\frac{S_{TE}^{\ast }}{K}){\displaystyle ⟨{\overline{\mathit{S}}}_T,{\mathbf{e}}_3⟩-({\alpha }_E+{\mu }_E){\overline{S}}_{TE}.}\end{array}$(52) ) is a decreasing function of $S_{TE}^{\ast }$ so that the unique solution to (51(51) $\begin{array}{r}0=r_E(1-\frac{S_{TE}^{\ast }}{K}){\displaystyle \frac{S_{TE}^{\ast }}{{\overline{S}}_{TE}}⟨{\overline{\mathit{S}}}_T,{\mathbf{e}}_3⟩-({\alpha }_E+{\mu }_E)S_{TE}^{\ast }}\end{array}$(51) ) is $S_{TE}^{\ast }={\overline{S}}_{TE}$. The equality (47(47) $\begin{array}{r}{\mathit{S}}_T^{\ast }+{\int }_0^{\mathrm{\infty }}{\mathit{i}}_T^{\ast }(a)\mathrm{d}a={\overline{\mathit{S}}}_T.\end{array}$(47) ) now follows from (50(50) $\begin{array}{r}{\int }_0^{\mathrm{\infty }}{\mathit{i}}_T^{\ast }(a)\mathrm{d}a+{\mathit{S}}_T^{\ast }={\alpha }_E{S}_{TE}^{\ast }{\mathbf{M}}^{-1}{\mathbf{e}}_1=\frac{S_{TE}^{\ast }}{{\overline{S}}_{TE}}{\overline{\mathit{S}}}_T.\end{array}$(50) ).

Thanks to Lemma 8.1 the system of equations (45(45) $\begin{array}{r}\{\begin{array}{l}{\displaystyle 0=r_E(1-\frac{S_{TE}^{\ast }}{K}){\displaystyle (⟨{\mathit{S}}_T^{\ast },{\mathbf{e}}_3⟩+{\int }_0^{\mathrm{\infty }}⟨{\mathit{i}}_T^{\ast }(a),{\mathbf{e}}_3⟩\mathrm{d}a)-({\alpha }_E+{\mu }_E)S_{TE}^{\ast }}}\\ {\displaystyle {\mathit{i}}_T^{\ast }(a)=\frac{1}{N_D^{\ast }}{\int }_0^{\mathrm{\infty }}⟨{\mathit{\beta }}_D(a),{\mathit{i}}_D^{\ast }(a)⟩\mathrm{d}a{e}^{-\mathbf{M}a}{\mathit{S}}_T^{\ast }}\\ {\displaystyle 0_{{\mathbb{R}}^3}={\alpha }_E{S}_{TE}^{\ast }{\mathbf{e}}_1-\frac{1}{N_D^{\ast }}{\int }_0^{\mathrm{\infty }}⟨{\mathit{\beta }}_D(a),{\mathit{i}}_D^{\ast }(a)⟩\mathrm{d}a{\mathit{S}}_T^{\ast }-\mathbf{M}{\mathit{S}}_T^{,}}\end{array}\end{array}$(45) ) is equivalent to (53) $\begin{array}{r}\{\begin{array}{l}{\displaystyle {\mathit{i}}_T^{\ast }(a)=\frac{1}{N_D^{\ast }}{\int }_0^{\mathrm{\infty }}⟨{\mathit{\beta }}_D(a),{\mathit{i}}_D^{\ast }(a)⟩\mathrm{d}a{e}^{-\mathbf{M}a}{\mathit{S}}_T^{\ast }}\\ {\displaystyle 0_{{\mathbb{R}}^3}={\alpha }_E{S}_{TE}{\mathbf{e}}_1-\frac{1}{N_D^{\ast }}{\int }_0^{\mathrm{\infty }}⟨{\mathit{\beta }}_D(a),{\mathit{i}}_D^{\ast }(a)⟩\mathrm{d}a{\mathit{S}}_T^{\ast }-\mathbf{M}{\mathit{S}}_T^{\ast }.}\end{array}\end{array}$(53) Note that for each nonnegative constant ${\eta }_0\ge 0$ the matrix $\mathbf{M}+{\eta }_0{\mathbf{I}}_3$ (where ${\mathbf{I}}_3$ is the identity operator on ${\mathbb{R}}^3$) is invertible with inverse $(\mathbf{M}+{\eta }_0{\mathbf{I}}_3{)}^{-1}$ satisfying (54) $\begin{array}{r}(\mathbf{M}+{\eta }_0{\mathbf{I}}_3{)}^{-1}{\mathbf{e}}_1=\sum _{k=1}^3\left(\prod _{j=1}^k\frac{{\alpha }_{j-1}}{{\mu }_j+{\alpha }_j+{\eta }_0}\right){\mathbf{e}}_k,\text{with }{\alpha }_0=1,{\alpha }_3=0,\end{array}$(54) and where the correspondence $L\leftarrow 1$, $N\leftarrow 2$, and $A\leftarrow 3$ is used. Therefore, using (53(53) $\begin{array}{r}\{\begin{array}{l}{\displaystyle {\mathit{i}}_T^{\ast }(a)=\frac{1}{N_D^{\ast }}{\int }_0^{\mathrm{\infty }}⟨{\mathit{\beta }}_D(a),{\mathit{i}}_D^{\ast }(a)⟩\mathrm{d}a{e}^{-\mathbf{M}a}{\mathit{S}}_T^{\ast }}\\ {\displaystyle 0_{{\mathbb{R}}^3}={\alpha }_E{S}_{TE}{\mathbf{e}}_1-\frac{1}{N_D^{\ast }}{\int }_0^{\mathrm{\infty }}⟨{\mathit{\beta }}_D(a),{\mathit{i}}_D^{\ast }(a)⟩\mathrm{d}a{\mathit{S}}_T^{\ast }-\mathbf{M}{\mathit{S}}_T^{\ast }.}\end{array}\end{array}$(53) ) together with (54(54) $\begin{array}{r}(\mathbf{M}+{\eta }_0{\mathbf{I}}_3{)}^{-1}{\mathbf{e}}_1=\sum _{k=1}^3\left(\prod _{j=1}^k\frac{{\alpha }_{j-1}}{{\mu }_j+{\alpha }_j+{\eta }_0}\right){\mathbf{e}}_k,\text{with }{\alpha }_0=1,{\alpha }_3=0,\end{array}$(54) ), for ${\eta }_0=\frac{1}{N_D^{\ast }}{\int }_0^{\mathrm{\infty }}⟨{\mathit{\beta }}_D(a),{\mathit{i}}_D^{\ast }(a)⟩\mathrm{d}a$, it follows that (55) $\begin{array}{r}{\mathit{S}}_T^{\ast }={\alpha }_E{S}_{TE}\sum _{k=1}^3{\mathbf{e}}_k\prod _{j=1}^k\frac{N_D^{\ast }{\alpha }_{j-1}}{N_D^{\ast }({\mu }_j+{\alpha }_j)+{\int }_0^{\mathrm{\infty }}⟨{\mathit{\beta }}_D(a),{\mathit{i}}_D^{\ast }(a)⟩\mathrm{d}a},\end{array}$(55) with ${\alpha }_0=1$ and ${\alpha }_3=0$. Thus, recalling the expression of ${\overline{\mathit{S}}}_T$ in (24(24) $\begin{array}{r}\{\begin{array}{l}{\overline{S}}_{TE}={\displaystyle K(1-\frac{1}{{\mathcal{R}}_{ST}^0}),}\\ {\displaystyle {\overline{\mathit{S}}}_T={\alpha }_E{\overline{S}}_{TE}\sum _{k=1}^3{\mathbf{e}}_k\prod _{j=1}^k\frac{{\alpha }_{j-1}}{{\mu }_j+{\alpha }_j},\text{with }{\alpha }_0=1,{\alpha }_3=0,}\end{array}\end{array}$(24) ) we obtain (56) $\begin{array}{r}{\mathit{S}}_T^{\ast }=\sum _{k=1}^3⟨{\overline{\mathit{S}}}_T,{\mathbf{e}}_k⟩{\mathbf{e}}_k\prod _{j=1}^k\frac{N_D^{\ast }({\mu }_j+{\alpha }_j)}{N_D^{\ast }({\mu }_j+{\alpha }_j)+{\int }_0^{\mathrm{\infty }}⟨{\mathit{\beta }}_D(a),{\mathit{i}}_D^{\ast }(a)⟩\mathrm{d}a}\end{array}$(56) and the system of equations (53(53) $\begin{array}{r}\{\begin{array}{l}{\displaystyle {\mathit{i}}_T^{\ast }(a)=\frac{1}{N_D^{\ast }}{\int }_0^{\mathrm{\infty }}⟨{\mathit{\beta }}_D(a),{\mathit{i}}_D^{\ast }(a)⟩\mathrm{d}a{e}^{-\mathbf{M}a}{\mathit{S}}_T^{\ast }}\\ {\displaystyle 0_{{\mathbb{R}}^3}={\alpha }_E{S}_{TE}{\mathbf{e}}_1-\frac{1}{N_D^{\ast }}{\int }_0^{\mathrm{\infty }}⟨{\mathit{\beta }}_D(a),{\mathit{i}}_D^{\ast }(a)⟩\mathrm{d}a{\mathit{S}}_T^{\ast }-\mathbf{M}{\mathit{S}}_T^{\ast }.}\end{array}\end{array}$(53) ) becomes equivalent to (57a) $\begin{array}{rl}& {\displaystyle {\mathit{S}}_T^{\ast }=\sum _{k=1}^3\left(\prod _{j=1}^k\frac{{\mu }_j+{\alpha }_j}{({\mu }_j+{\alpha }_j)+W_D^{\ast }}\right)⟨{\overline{\mathit{S}}}_T,{\mathbf{e}}_k⟩{\mathbf{e}}_k}\end{array}$(57a) (57b) $\begin{array}{rl}& {\displaystyle {\mathit{i}}_T^{\ast }(a)=W_D^{\diamond }{e}^{-\mathbf{M}a}{\mathit{S}}_T^{\ast }}\end{array}$(57b) (57c) $\begin{array}{rl}& {\displaystyle W_D^{\diamond }={\int }_0^{\mathrm{\infty }}⟨{\mathit{\beta }}_D(a),\frac{{\mathit{i}}_D^{\ast }(a)}{N_D^{\ast }}⟩\mathrm{d}a.}\end{array}$(57c) From the above comments, it follows that ${\mathit{S}}_T^{\ast }$ and ${\mathit{i}}_T^{\ast }$ are entirely determined by the variables $W_D^{\diamond }$ and $\frac{{\mathit{i}}_D^{\ast }}{N_D^{\ast }}$.

On the stationary states of dogs: In order to obtain an equivalent formulation to (46(46) $\begin{array}{r}\{\begin{array}{l}{\displaystyle 0={\theta }_D-\frac{S_D^{\ast }}{N_D^{\ast }}{\int }_0^{\mathrm{\infty }}{\beta }_T(a)⟨\mathbf{1},{\mathit{i}}_T^{\ast }(a)⟩\mathrm{d}a-{\mu }_D{S}_D^{\ast }}\\ {\displaystyle {\mathit{i}}_D^{\ast }(a)=\frac{S_D^{\ast }+ϵR_D^{\ast }}{N_D^{\ast }}{\int }_0^{\mathrm{\infty }}{\beta }_T(a)⟨\mathbf{1},{\mathit{i}}_T^{\ast }(a)⟩\mathrm{d}a{\pi }_D(a){\mathbf{v}}_1}\\ {\displaystyle +{\int }_0^{\mathrm{\infty }}{\nu }_D(\tau )⟨{\mathbf{v}}_1,{\mathit{i}}_D^{\ast }(a)⟩\mathrm{d}a{\pi }_C(a){\mathbf{v}}_2}\\ {\displaystyle 0={\int }_0^{\mathrm{\infty }}\gamma (a)⟨{\mathbf{v}}_1,{\mathit{i}}_D^{\ast }(a)⟩\mathrm{d}a-ϵ\frac{R_D^{\ast }}{N_D^{\ast }}{\int }_0^{\mathrm{\infty }}{\beta }_T(a)⟨\mathbf{1},{\mathit{i}}_T^{\ast }(a)⟩\mathrm{d}a-{\mu }_D{R}_D^{\ast }}\\ {\displaystyle 0={\theta }_D-{\mu }_D{N}_D^{\ast }-{\int }_0^{\mathrm{\infty }}{\delta }_D(a)⟨{\mathbf{v}}_1,{\mathit{i}}_D^{\ast }(a)⟩\mathrm{d}a-{\int }_0^{\mathrm{\infty }}{\delta }_C(a)⟨{\mathbf{v}}_2,{\mathit{i}}_D^{\ast }(a)⟩\mathrm{d}a.}\end{array}\end{array}$(46) ) we introduce the new variables (58) $\begin{array}{r}{\mathit{i}}_D^{\diamond }(a):=\frac{{\mathit{i}}_D^{\ast }(a)}{N_D^{\ast }},S_D^{\diamond }:=\frac{S_D^{\ast }}{N_D^{\ast }}\mathrm{and}{R}_D^{\diamond }:=\frac{R_D^{\ast }}{N_D^{\ast }}.\end{array}$(58) Next, we set (59) $\begin{array}{r}{P}_T^{\ast }:={\int }_0^{\mathrm{\infty }}{\beta }_T(a)⟨\mathbf{1},{\mathit{i}}_T^{\ast }(a)⟩\mathrm{d}a,\end{array}$(59) and (60) $\begin{array}{r}\{\begin{array}{l}{\displaystyle P_D^{\diamond }:={\int }_0^{\mathrm{\infty }}\gamma (a)⟨{\mathbf{v}}_1,{\mathit{i}}_D^{\diamond }(a)⟩\mathrm{d}a}\\ {\displaystyle Q_D^{\diamond }:={\int }_0^{\mathrm{\infty }}{\nu }_D(a)⟨{\mathbf{v}}_1,{\mathit{i}}_D^{\diamond }(a)⟩\mathrm{d}a,}\\ {\displaystyle U_D^{\diamond }:={\int }_0^{\mathrm{\infty }}{\delta }_D(a)⟨{\mathbf{v}}_1,{\mathit{i}}_D^{\ast }(a)⟩\mathrm{d}a}\\ {\displaystyle U_C^{\diamond }:={\int }_0^{\mathrm{\infty }}{\delta }_C(a)⟨{\mathbf{v}}_2,{\mathit{i}}_D^{\ast }(a)⟩\mathrm{d}a}\end{array}\end{array}$(60) With these new variables, system (46(46) $\begin{array}{r}\{\begin{array}{l}{\displaystyle 0={\theta }_D-\frac{S_D^{\ast }}{N_D^{\ast }}{\int }_0^{\mathrm{\infty }}{\beta }_T(a)⟨\mathbf{1},{\mathit{i}}_T^{\ast }(a)⟩\mathrm{d}a-{\mu }_D{S}_D^{\ast }}\\ {\displaystyle {\mathit{i}}_D^{\ast }(a)=\frac{S_D^{\ast }+ϵR_D^{\ast }}{N_D^{\ast }}{\int }_0^{\mathrm{\infty }}{\beta }_T(a)⟨\mathbf{1},{\mathit{i}}_T^{\ast }(a)⟩\mathrm{d}a{\pi }_D(a){\mathbf{v}}_1}\\ {\displaystyle +{\int }_0^{\mathrm{\infty }}{\nu }_D(\tau )⟨{\mathbf{v}}_1,{\mathit{i}}_D^{\ast }(a)⟩\mathrm{d}a{\pi }_C(a){\mathbf{v}}_2}\\ {\displaystyle 0={\int }_0^{\mathrm{\infty }}\gamma (a)⟨{\mathbf{v}}_1,{\mathit{i}}_D^{\ast }(a)⟩\mathrm{d}a-ϵ\frac{R_D^{\ast }}{N_D^{\ast }}{\int }_0^{\mathrm{\infty }}{\beta }_T(a)⟨\mathbf{1},{\mathit{i}}_T^{\ast }(a)⟩\mathrm{d}a-{\mu }_D{R}_D^{\ast }}\\ {\displaystyle 0={\theta }_D-{\mu }_D{N}_D^{\ast }-{\int }_0^{\mathrm{\infty }}{\delta }_D(a)⟨{\mathbf{v}}_1,{\mathit{i}}_D^{\ast }(a)⟩\mathrm{d}a-{\int }_0^{\mathrm{\infty }}{\delta }_C(a)⟨{\mathbf{v}}_2,{\mathit{i}}_D^{\ast }(a)⟩\mathrm{d}a.}\end{array}\end{array}$(46) ) becomes (61) $\begin{array}{r}\{\begin{array}{l}{\displaystyle 0={\theta }_D-S_D^{\diamond }{P}_T^{\ast }-{\mu }_D{N}_D^{\ast }{S}_D^{\diamond }}\\ {\displaystyle N_D^{\ast }{\mathit{i}}_D^{\diamond }(a)=(S_D^{\diamond }+ϵR_D^{\diamond })P_T^{\ast }{\pi }_D(a){\mathbf{v}}_1+N_D^{\ast }{Q}_D^{\diamond }{\pi }_C(a){\mathbf{v}}_2}\\ {\displaystyle 0=N_D^{\ast }{P}_D^{\diamond }-ϵR_D^{\diamond }{P}_T^{\ast }-{\mu }_D{N}_D^{\ast }{R}_D^{\diamond }}\\ 0={\theta }_D-{\mu }_D{N}_D^{\ast }-N_D^{\ast }{U}_D^{\diamond }-N_D^{\ast }{U}_C^{\diamond }\end{array}\end{array}$(61) which is equivalent to (62) $\begin{array}{r}\{\begin{array}{l}{\displaystyle N_D^{\ast }=\frac{{\theta }_D}{{\mu }_D+U_D^{\diamond }+U_C^{\diamond }}}\\ {\displaystyle S_D^{\diamond }=\frac{{\theta }_D}{{\mu }_D{N}_D^{\ast }+P_T^{\ast }}}\\ {\displaystyle R_D^{\diamond }=\frac{N_D^{\ast }{P}_D^{\diamond }}{ϵP_T^{\ast }+{\mu }_D{N}_D^{\ast }}}\\ {\displaystyle {\mathit{i}}_D^{\diamond }(a)=\frac{S_D^{\diamond }+ϵR_D^{\diamond }}{N_D^{\ast }}{P}_T^{\ast }{\pi }_D(a){\mathbf{v}}_1+Q_D^{\diamond }{\pi }_C(a){\mathbf{v}}_2.}\end{array}\end{array}$(62) Next, we show the relationship between $P_T^{\ast }$ and $W_D^{\diamond }$ defined respectively in (57c(57c) $\begin{array}{rl}& {\displaystyle W_D^{\diamond }={\int }_0^{\mathrm{\infty }}⟨{\mathit{\beta }}_D(a),\frac{{\mathit{i}}_D^{\ast }(a)}{N_D^{\ast }}⟩\mathrm{d}a.}\end{array}$(57c) ) and (59(59) $\begin{array}{r}{P}_T^{\ast }:={\int }_0^{\mathrm{\infty }}{\beta }_T(a)⟨\mathbf{1},{\mathit{i}}_T^{\ast }(a)⟩\mathrm{d}a,\end{array}$(59) ). To do so, we integrate (57b(57b) $\begin{array}{rl}& {\displaystyle {\mathit{i}}_T^{\ast }(a)=W_D^{\diamond }{e}^{-\mathbf{M}a}{\mathit{S}}_T^{\ast }}\end{array}$(57b) ) to obtain (63) $\begin{array}{r}{P}_T^{\ast }={\displaystyle W_D^{\diamond }{\int }_0^{\mathrm{\infty }}{\beta }_T(a)⟨\mathbf{1},e^{-\mathbf{M}a}{\mathit{S}}_T^{\ast }⟩\mathrm{d}a}\end{array}$(63) and by using the expression of ${\mathit{S}}_T^{\ast }$ defined in (57a(57a) $\begin{array}{rl}& {\displaystyle {\mathit{S}}_T^{\ast }=\sum _{k=1}^3\left(\prod _{j=1}^k\frac{{\mu }_j+{\alpha }_j}{({\mu }_j+{\alpha }_j)+W_D^{\ast }}\right)⟨{\overline{\mathit{S}}}_T,{\mathbf{e}}_k⟩{\mathbf{e}}_k}\end{array}$(57a) ) we obtain the following more explicit formula (64) $\begin{array}{r}{P}_T^{\ast }=W_D^{\diamond }\sum _{k=1}^3\left(\prod _{j=1}^k\frac{{\mu }_j+{\alpha }_j}{({\mu }_j+{\alpha }_j)+W_D^{\diamond }}\right)\frac{{\overline{S}}_k}{\sum _{j\in \mathcal{S}\setminus \{E\}}{\overline{S}}_j}{\mathcal{T}}_{0k}^{T\to D}{\overline{S}}_D,\end{array}$(64) with ${\mathcal{T}}_{0k}^{T\to D}$ the k's stage vectorial capacity defined in (17(17) $\begin{array}{r}{\mathcal{T}}_{0k}^{T\to D}=\underset{\begin{array}{c}\text{Ticks/Dogs ratio}\end{array}}{\underbrace{\frac{\sum _{j\in \mathcal{S}\setminus \{E\}}{\overline{S}}_j}{{\overline{S}}_D}}}\times \underset{\begin{array}{c}\text{Transmission from ticks at stage }k\text{ to Dogs}\end{array}}{\underbrace{{\int }_0^{\mathrm{\infty }}{\beta }_T(a)⟨\mathbf{1},e^{-\mathbf{M}a}{\mathbf{e}}_k⟩\mathrm{d}a,}}\end{array}$(17) ).

In the following, we determine the equations for the variables $P_D^{\diamond }$, $Q_D^{\diamond }$, $U_D^{\diamond }$, $U_C^{\diamond }$ and $W_C^{\diamond }$ defined in (60(60) $\begin{array}{r}\{\begin{array}{l}{\displaystyle P_D^{\diamond }:={\int }_0^{\mathrm{\infty }}\gamma (a)⟨{\mathbf{v}}_1,{\mathit{i}}_D^{\diamond }(a)⟩\mathrm{d}a}\\ {\displaystyle Q_D^{\diamond }:={\int }_0^{\mathrm{\infty }}{\nu }_D(a)⟨{\mathbf{v}}_1,{\mathit{i}}_D^{\diamond }(a)⟩\mathrm{d}a,}\\ {\displaystyle U_D^{\diamond }:={\int }_0^{\mathrm{\infty }}{\delta }_D(a)⟨{\mathbf{v}}_1,{\mathit{i}}_D^{\ast }(a)⟩\mathrm{d}a}\\ {\displaystyle U_C^{\diamond }:={\int }_0^{\mathrm{\infty }}{\delta }_C(a)⟨{\mathbf{v}}_2,{\mathit{i}}_D^{\ast }(a)⟩\mathrm{d}a}\end{array}\end{array}$(60) ). To this end, recalling that ${\mathit{\beta }}_D(a)={\beta }_D^1(a){\mathbf{v}}_1+{\beta }_D^2(a){\mathbf{v}}_1$, we use successively the equation of ${\mathit{i}}_D^{\diamond }(a)$ in (62(62) $\begin{array}{r}\{\begin{array}{l}{\displaystyle N_D^{\ast }=\frac{{\theta }_D}{{\mu }_D+U_D^{\diamond }+U_C^{\diamond }}}\\ {\displaystyle S_D^{\diamond }=\frac{{\theta }_D}{{\mu }_D{N}_D^{\ast }+P_T^{\ast }}}\\ {\displaystyle R_D^{\diamond }=\frac{N_D^{\ast }{P}_D^{\diamond }}{ϵP_T^{\ast }+{\mu }_D{N}_D^{\ast }}}\\ {\displaystyle {\mathit{i}}_D^{\diamond }(a)=\frac{S_D^{\diamond }+ϵR_D^{\diamond }}{N_D^{\ast }}{P}_T^{\ast }{\pi }_D(a){\mathbf{v}}_1+Q_D^{\diamond }{\pi }_C(a){\mathbf{v}}_2.}\end{array}\end{array}$(62) ) to obtain (65a) $\begin{array}{rl}{\displaystyle W_D^{\diamond }}& =\frac{S_D^{\diamond }+ϵR_D^{\diamond }}{N_D^{\ast }}{P}_T^{\ast }{\int }_0^{\mathrm{\infty }}{\beta }_D^1(a){\pi }_D(a)\mathrm{d}a\end{array}$(65a) (65b) $\begin{array}{rl}& {\displaystyle +Q_D^{\diamond }{\int }_0^{\mathrm{\infty }}{\beta }_D^2(a){\pi }_C(a)\mathrm{d}a}\end{array}$(65b) (65c) $\begin{array}{rl}{\displaystyle P_D^{\diamond }}& =\frac{S_D^{\diamond }+ϵR_D^{\diamond }}{N_D^{\ast }}{P}_T^{\ast }{\int }_0^{\mathrm{\infty }}\gamma (a){\pi }_D(a)\mathrm{d}a\end{array}$(65c) (65d) $\begin{array}{rl}{\displaystyle Q_D^{\diamond }}& =\frac{S_D^{\diamond }+ϵR_D^{\diamond }}{N_D^{\ast }}{P}_T^{\ast }{\int }_0^{\mathrm{\infty }}{\nu }_D(a){\pi }_D(a)\mathrm{d}a\end{array}$(65d) (65e) $\begin{array}{rl}{\displaystyle U_D^{\diamond }}& =\frac{S_D^{\diamond }+ϵR_D^{\diamond }}{N_D^{\ast }}{P}_T^{\ast }{\int }_0^{\mathrm{\infty }}{\delta }_D(a){\pi }_D(a)\mathrm{d}a\end{array}$(65e) (65f) $\begin{array}{rl}{\displaystyle U_C^{\diamond }}& =\frac{S_D^{\diamond }+ϵR_D^{\diamond }}{N_D^{\ast }}{P}_T^{\ast }{\int }_0^{\mathrm{\infty }}{\delta }_C(a){\pi }_C(a)\mathrm{d}a.\end{array}$(65f)

Note that, plugging (65d(65d) $\begin{array}{rl}{\displaystyle Q_D^{\diamond }}& =\frac{S_D^{\diamond }+ϵR_D^{\diamond }}{N_D^{\ast }}{P}_T^{\ast }{\int }_0^{\mathrm{\infty }}{\nu }_D(a){\pi }_D(a)\mathrm{d}a\end{array}$(65d) ) into (65a(65a) $\begin{array}{rl}{\displaystyle W_D^{\diamond }}& =\frac{S_D^{\diamond }+ϵR_D^{\diamond }}{N_D^{\ast }}{P}_T^{\ast }{\int }_0^{\mathrm{\infty }}{\beta }_D^1(a){\pi }_D(a)\mathrm{d}a\end{array}$(65a) ) gives the following equality (66) $\begin{array}{r}\{\begin{array}{l}{\displaystyle W_D^{\diamond }=\frac{S_D^{\diamond }+ϵR_D^{\diamond }}{N_D^{\ast }}{P}_T^{\ast }[{\int }_0^{\mathrm{\infty }}{\beta }_D^1(a){\pi }_D(a)\mathrm{d}a}\\ {\displaystyle +{\int }_0^{\mathrm{\infty }}{\nu }_D(a){\pi }_D(a)\mathrm{d}a{\int }_0^{\mathrm{\infty }}{\beta }_D^2(a){\pi }_C(a)\mathrm{d}a].}\end{array}\end{array}$(66) Let us now define the function (67) $\begin{array}{r}f(W_D^{\diamond }):=\sum _{k=1}^3\left(\prod _{j=1}^k\frac{{\mu }_j+{\alpha }_j}{({\mu }_j+{\alpha }_j)+W_D^{\diamond }}\right)\frac{{\overline{S}}_k}{\sum _{j\in \mathcal{S}\setminus \{E\}}{\overline{S}}_j}{\mathcal{T}}_{0k}^{T\to D}{\overline{S}}_D\end{array}$(67) and observe that the vectorial capacity ${\mathcal{T}}_0^{T\to D}$ is given by (68) $\begin{array}{r}{\mathcal{T}}_0^{T\to D}=\frac{f(0)}{{\overline{S}}_D}\end{array}$(68) so that (69) $\begin{array}{r}\{\begin{array}{l}{\displaystyle {\mathcal{T}}_0=\frac{f(0)}{{\overline{S}}_D}[{\int }_0^{\mathrm{\infty }}{\beta }_D^1(a){\pi }_D(a)\mathrm{d}a}\\ {\displaystyle +{\int }_0^{\mathrm{\infty }}{\nu }_D(a){\pi }_D(a)\mathrm{d}a{\int }_0^{\mathrm{\infty }}{\beta }_D^2(a){\pi }_C(a)\mathrm{d}a].}\end{array}\end{array}$(69) In the following, we rewrite our system of equations by using the variable (70) $\begin{array}{r}K:=\frac{S_D^{\diamond }+ϵR_D^{\diamond }}{N_D^{\ast }}\frac{P_T^{\ast }}{{\mathcal{T}}_0^{T\to D}}.\end{array}$(70) To do so, we first observe that (66(66) $\begin{array}{r}\{\begin{array}{l}{\displaystyle W_D^{\diamond }=\frac{S_D^{\diamond }+ϵR_D^{\diamond }}{N_D^{\ast }}{P}_T^{\ast }[{\int }_0^{\mathrm{\infty }}{\beta }_D^1(a){\pi }_D(a)\mathrm{d}a}\\ {\displaystyle +{\int }_0^{\mathrm{\infty }}{\nu }_D(a){\pi }_D(a)\mathrm{d}a{\int }_0^{\mathrm{\infty }}{\beta }_D^2(a){\pi }_C(a)\mathrm{d}a].}\end{array}\end{array}$(66) ) and (64(64) $\begin{array}{r}{P}_T^{\ast }=W_D^{\diamond }\sum _{k=1}^3\left(\prod _{j=1}^k\frac{{\mu }_j+{\alpha }_j}{({\mu }_j+{\alpha }_j)+W_D^{\diamond }}\right)\frac{{\overline{S}}_k}{\sum _{j\in \mathcal{S}\setminus \{E\}}{\overline{S}}_j}{\mathcal{T}}_{0k}^{T\to D}{\overline{S}}_D,\end{array}$(64) ) take respectively the following form (71) $\begin{array}{r}{W}_D^{\diamond }=K{\mathcal{T}}_0\mathrm{and}{P}_T^{\ast }=W_D^{\diamond }f(W_D^{\diamond })=K{\mathcal{T}}_{0}f(K{\mathcal{T}}_0)\end{array}$(71) while (65a(65a) $\begin{array}{rl}{\displaystyle W_D^{\diamond }}& =\frac{S_D^{\diamond }+ϵR_D^{\diamond }}{N_D^{\ast }}{P}_T^{\ast }{\int }_0^{\mathrm{\infty }}{\beta }_D^1(a){\pi }_D(a)\mathrm{d}a\end{array}$(65a) ) is now given by (72a) $\begin{array}{rl}{\displaystyle W_D^{\diamond }}& =\frac{S_D^{\diamond }+ϵR_D^{\diamond }}{N_D^{\ast }}{P}_T^{\ast }{\int }_0^{\mathrm{\infty }}{\beta }_D^1(a){\pi }_D(a)\mathrm{d}a\end{array}$(72a) (72b) $\begin{array}{rl}& {\displaystyle +Q_D^{\diamond }{\int }_0^{\mathrm{\infty }}{\beta }_D^2(a){\pi }_C(a)\mathrm{d}a}\end{array}$(72b) (72c) $\begin{array}{rl}{\displaystyle P_D^{\diamond }}& =K{\mathcal{T}}_0^{T\to D}{\int }_0^{\mathrm{\infty }}\gamma (a){\pi }_D(a)\mathrm{d}a\end{array}$(72c) (72d) $\begin{array}{rl}{\displaystyle Q_D^{\diamond }}& =K{\mathcal{T}}_0^{T\to D}{\int }_0^{\mathrm{\infty }}{\nu }_D(a){\pi }_D(a)\mathrm{d}a\end{array}$(72d) (72e) $\begin{array}{rl}{\displaystyle U_D^{\diamond }}& =K{\mathcal{T}}_0^{T\to D}{\int }_0^{\mathrm{\infty }}{\delta }_D(a){\pi }_D(a)\mathrm{d}a\end{array}$(72e) (72f) $\begin{array}{rl}{\displaystyle U_C^{\diamond }}& =K{\mathcal{T}}_0^{T\to D}{\int }_0^{\mathrm{\infty }}{\delta }_C(a){\pi }_C(a)\mathrm{d}a.\end{array}$(72f) Hence, using the first equation of (62(62) $\begin{array}{r}\{\begin{array}{l}{\displaystyle N_D^{\ast }=\frac{{\theta }_D}{{\mu }_D+U_D^{\diamond }+U_C^{\diamond }}}\\ {\displaystyle S_D^{\diamond }=\frac{{\theta }_D}{{\mu }_D{N}_D^{\ast }+P_T^{\ast }}}\\ {\displaystyle R_D^{\diamond }=\frac{N_D^{\ast }{P}_D^{\diamond }}{ϵP_T^{\ast }+{\mu }_D{N}_D^{\ast }}}\\ {\displaystyle {\mathit{i}}_D^{\diamond }(a)=\frac{S_D^{\diamond }+ϵR_D^{\diamond }}{N_D^{\ast }}{P}_T^{\ast }{\pi }_D(a){\mathbf{v}}_1+Q_D^{\diamond }{\pi }_C(a){\mathbf{v}}_2.}\end{array}\end{array}$(62) ) together with (72e(72e) $\begin{array}{rl}{\displaystyle U_D^{\diamond }}& =K{\mathcal{T}}_0^{T\to D}{\int }_0^{\mathrm{\infty }}{\delta }_D(a){\pi }_D(a)\mathrm{d}a\end{array}$(72e) ), and (72f(72f) $\begin{array}{rl}{\displaystyle U_C^{\diamond }}& =K{\mathcal{T}}_0^{T\to D}{\int }_0^{\mathrm{\infty }}{\delta }_C(a){\pi }_C(a)\mathrm{d}a.\end{array}$(72f) ) we obtain the following necessary conditions (73) $\begin{array}{r}{N}_D^{\ast }(K):=\frac{{\theta }_D}{{\mu }_D+K{\mathcal{T}}_0^{T\to D}({\int }_0^{\mathrm{\infty }}{\delta }_D(a){\pi }_D(a)\mathrm{d}a+{\int }_0^{\mathrm{\infty }}{\delta }_C(a){\pi }_C(a)\mathrm{d}a)}.\end{array}$(73) Furthermore, the second, third and fourth equation of (62(62) $\begin{array}{r}\{\begin{array}{l}{\displaystyle N_D^{\ast }=\frac{{\theta }_D}{{\mu }_D+U_D^{\diamond }+U_C^{\diamond }}}\\ {\displaystyle S_D^{\diamond }=\frac{{\theta }_D}{{\mu }_D{N}_D^{\ast }+P_T^{\ast }}}\\ {\displaystyle R_D^{\diamond }=\frac{N_D^{\ast }{P}_D^{\diamond }}{ϵP_T^{\ast }+{\mu }_D{N}_D^{\ast }}}\\ {\displaystyle {\mathit{i}}_D^{\diamond }(a)=\frac{S_D^{\diamond }+ϵR_D^{\diamond }}{N_D^{\ast }}{P}_T^{\ast }{\pi }_D(a){\mathbf{v}}_1+Q_D^{\diamond }{\pi }_C(a){\mathbf{v}}_2.}\end{array}\end{array}$(62) ) together with (71(71) $\begin{array}{r}{W}_D^{\diamond }=K{\mathcal{T}}_0\mathrm{and}{P}_T^{\ast }=W_D^{\diamond }f(W_D^{\diamond })=K{\mathcal{T}}_{0}f(K{\mathcal{T}}_0)\end{array}$(71) ), (72c(72c) $\begin{array}{rl}{\displaystyle P_D^{\diamond }}& =K{\mathcal{T}}_0^{T\to D}{\int }_0^{\mathrm{\infty }}\gamma (a){\pi }_D(a)\mathrm{d}a\end{array}$(72c) ) and (72d(72d) $\begin{array}{rl}{\displaystyle Q_D^{\diamond }}& =K{\mathcal{T}}_0^{T\to D}{\int }_0^{\mathrm{\infty }}{\nu }_D(a){\pi }_D(a)\mathrm{d}a\end{array}$(72d) ) provide (74a) $\begin{array}{rl}{\displaystyle R_D^{\diamond }({\mathcal{T}}_0,K)}& :={\mathcal{T}}_0^{T\to D}\frac{K{N}_D^{\ast }(K)}{ϵK{\mathcal{T}}_{0}f(K{\mathcal{T}}_0)+{\mu }_D{N}_D^{\ast }(K)}{\int }_0^{\mathrm{\infty }}\gamma (a){\pi }_D(a)\mathrm{d}a\end{array}$(74a) (74b) $\begin{array}{rl}{\displaystyle S_D^{\diamond }({\mathcal{T}}_0,K)}& :=\frac{{\theta }_D}{{\mu }_D{N}_D^{\ast }(K)+K{\mathcal{T}}_{0}f(K{\mathcal{T}}_0)}\end{array}$(74b) (74c) $\begin{array}{rl}{\displaystyle {\mathit{i}}_D^{\diamond }(a)}& :=K{\mathcal{T}}_0^{T\to D}{\pi }_D(a){\mathbf{v}}_1+K{\mathcal{T}}_0^{T\to D}{\int }_0^{\mathrm{\infty }}{\nu }_D(a){\pi }_D(a)\mathrm{d}a{\pi }_C(a){\mathbf{v}}_2.\end{array}$(74c) Moreover, from the equality $P_T^{\ast }=K{\mathcal{T}}_{0}f(K{\mathcal{T}}_0)$ and (70(70) $\begin{array}{r}K:=\frac{S_D^{\diamond }+ϵR_D^{\diamond }}{N_D^{\ast }}\frac{P_T^{\ast }}{{\mathcal{T}}_0^{T\to D}}.\end{array}$(70) ) it comes that for K>0 we must have the following equality (75) $\begin{array}{r}1=\frac{S_D^{\diamond }({\mathcal{T}}_0,K)+ϵR_D^{\diamond }({\mathcal{T}}_0,K)}{N_D^{\ast }(K)}\frac{f(K{\mathcal{T}}_0)}{{\mathcal{T}}_0^{T\to D}}{\mathcal{T}}_0.\end{array}$(75) Recalling that we have the necessary condition $S_D^{\diamond }+R_D^{\diamond }+{\int }_0^{\mathrm{\infty }}⟨\mathbf{1},{\mathit{i}}_D^{\diamond }(a)⟩\mathrm{d}a=1$, we consider the following equation $\begin{array}{r}\mathrm{\Delta }({\mathcal{T}}_0,K)=1\end{array}$with (76) $\begin{array}{r}\mathrm{\Delta }({\mathcal{T}}_0,K):=\frac{1-(1-ϵ)R_D^{\diamond }({\mathcal{T}}_0,K)-I_D^{\diamond }({\mathcal{T}}_0,K)}{N_D^{\ast }(K)}\frac{f(K{\mathcal{T}}_0)}{{\mathcal{T}}_0^{T\to D}}{\mathcal{T}}_0\end{array}$(76) and (77) $\begin{array}{r}{I}_D^{\diamond }({\mathcal{T}}_0,K):=K{\mathcal{T}}_0^{T\to D}{\int }_0^{\mathrm{\infty }}{\pi }_D(a)\mathrm{d}a+K{\mathcal{T}}_0^{T\to D}{\int }_0^{\mathrm{\infty }}{\nu }_D(a){\pi }_D(a)\mathrm{d}a{\int }_0^{\mathrm{\infty }}{\pi }_C(a)\mathrm{d}a.\end{array}$(77)

Lemma 8.2

There exists an endemic equilibrium to (7(7) $\begin{array}{r}\{\begin{array}{l}{\displaystyle \frac{\mathrm{d}{S}_D(t)}{\mathrm{d}t}={\theta }_D-\frac{S_D(t)}{N_D(t)}{\int }_0^{\mathrm{\infty }}{\beta }_T(a)⟨\mathbf{1},{\mathit{i}}_T(t,a)⟩\mathrm{d}a-{\mu }_D{S}_D(t),}\\ {\displaystyle \frac{\mathrm{\partial }{\mathit{i}}_D(t,a)}{\mathrm{\partial }t}+\frac{\mathrm{\partial }{\mathit{i}}_D(t,a)}{\mathrm{\partial }a}=-({\mu }_D+\mathit{\theta }(a)){\mathit{i}}_D(t,a),}\\ {\mathit{i}}_D(t,0)={\displaystyle \frac{S_D(t)+ϵR_D(t)}{N_D(t)}{\int }_0^{\mathrm{\infty }}{\beta }_T(a)⟨\mathbf{1},{\mathit{i}}_T(t,a)⟩\mathrm{d}a{\mathbf{v}}_1}\\ {\displaystyle +{\int }_0^{\mathrm{\infty }}{\nu }_D(a)⟨{\mathbf{v}}_1,{\mathit{i}}_D(t,a)⟩\mathrm{d}a{\mathbf{v}}_2,}\\ {\displaystyle \frac{\mathrm{d}{R}_D(t)}{\mathrm{d}t}={\int }_0^{\mathrm{\infty }}\gamma (a)⟨{\mathbf{v}}_1,{\mathit{i}}_D(t,a)⟩\mathrm{d}a-ϵ\frac{R_D(t)}{N_D(t)}{\int }_0^{\mathrm{\infty }}{\beta }_T(a)⟨\mathbf{1},{\mathit{i}}_T(t,a)⟩\mathrm{d}a-{\mu }_D{R}_D(t),}\end{array}\end{array}$(7) )–(8(8) $\begin{array}{r}\{\begin{array}{l}{\displaystyle \frac{\mathrm{d}{S}_{TE}(t)}{\mathrm{d}t}=r_E(1-\frac{S_{TE}(t)}{K}){\displaystyle (⟨{\mathit{S}}_T(t),{\mathbf{e}}_3⟩+{\int }_0^{\mathrm{\infty }}⟨{\mathit{i}}_T(t,a),{\mathbf{e}}_3⟩\mathrm{d}a)-({\alpha }_E+{\mu }_E)S_{TE}(t),}}\\ {\displaystyle \frac{\mathrm{d}{\mathit{S}}_T(t)}{\mathrm{d}t}={\alpha }_E{S}_{TE}(t){\mathbf{e}}_1-\mathbf{M}{\mathit{S}}_T(t)-\frac{{\mathit{S}}_T(t)}{N_D(t)}{\int }_0^{\mathrm{\infty }}⟨{\mathit{\beta }}_D(a),{\mathit{i}}_D(t,a)⟩\mathrm{d}a,}\\ {\displaystyle \frac{\mathrm{\partial }{\mathit{i}}_T(t,a)}{\mathrm{\partial }t}+\frac{\mathrm{\partial }{\mathit{i}}_T(t,a)}{\mathrm{\partial }a}=-\mathbf{M}{\mathit{i}}_T(t,a),}\\ {\mathit{i}}_T(t,0)={\displaystyle \frac{{\mathit{S}}_T(t)}{N_D(t)}{\int }_0^{\mathrm{\infty }}⟨{\mathit{\beta }}_D(a),{\mathit{i}}_D(t,a)⟩\mathrm{d}a,}\end{array}\end{array}$(8) ) if and only if there exists K>0 such that $\mathrm{\Delta }({\mathcal{T}}_0,K)=1$.

Proof.

The necessity of the lemma follows from the preceding arguments. Next, we prove the sufficiency. To this end, assume that there exists K>0 such that $\mathrm{\Delta }({\mathcal{T}}_0,K)=1$. Let $N_D^{\ast }$ be given by the right-hand side of (73(73) $\begin{array}{r}{N}_D^{\ast }(K):=\frac{{\theta }_D}{{\mu }_D+K{\mathcal{T}}_0^{T\to D}({\int }_0^{\mathrm{\infty }}{\delta }_D(a){\pi }_D(a)\mathrm{d}a+{\int }_0^{\mathrm{\infty }}{\delta }_C(a){\pi }_C(a)\mathrm{d}a)}.\end{array}$(73) ). Let $R_D^{\diamond }$ and ${\mathit{i}}_D^{\diamond }$ be given respectively by the right-hand side of (74a(74a) $\begin{array}{rl}{\displaystyle R_D^{\diamond }({\mathcal{T}}_0,K)}& :={\mathcal{T}}_0^{T\to D}\frac{K{N}_D^{\ast }(K)}{ϵK{\mathcal{T}}_{0}f(K{\mathcal{T}}_0)+{\mu }_D{N}_D^{\ast }(K)}{\int }_0^{\mathrm{\infty }}\gamma (a){\pi }_D(a)\mathrm{d}a\end{array}$(74a) ) and (74c(74c) $\begin{array}{rl}{\displaystyle {\mathit{i}}_D^{\diamond }(a)}& :=K{\mathcal{T}}_0^{T\to D}{\pi }_D(a){\mathbf{v}}_1+K{\mathcal{T}}_0^{T\to D}{\int }_0^{\mathrm{\infty }}{\nu }_D(a){\pi }_D(a)\mathrm{d}a{\pi }_C(a){\mathbf{v}}_2.\end{array}$(74c) ). Observe that with these definitions we have (78) $\begin{array}{r}{I}_D^{\diamond }({\mathcal{T}}_0,K)={\int }_0^{\mathrm{\infty }}⟨\mathbf{1},{\mathit{i}}_D^{\diamond }⟩\mathrm{d}a\end{array}$(78) so that setting (79) $\begin{array}{r}{S}_D^{\diamond }:=1-R_D^{\diamond }-I_D^{\diamond }({\mathcal{T}}_0,K)\end{array}$(79) it comes from (76(76) $\begin{array}{r}\mathrm{\Delta }({\mathcal{T}}_0,K):=\frac{1-(1-ϵ)R_D^{\diamond }({\mathcal{T}}_0,K)-I_D^{\diamond }({\mathcal{T}}_0,K)}{N_D^{\ast }(K)}\frac{f(K{\mathcal{T}}_0)}{{\mathcal{T}}_0^{T\to D}}{\mathcal{T}}_0\end{array}$(76) ) (80) $\begin{array}{r}1=\frac{S_D^{\diamond }+ϵR_D^{\diamond }}{N_D^{\ast }}\frac{f(K{\mathcal{T}}_0)}{{\mathcal{T}}_0^{T\to D}}{\mathcal{T}}_0.\end{array}$(80) Note that setting $W_D^{\diamond }=K{\mathcal{T}}_0$, $P_T^{\ast }=W_D^{\diamond }f(W_D^{\diamond })$, with f defined in (67(67) $\begin{array}{r}f(W_D^{\diamond }):=\sum _{k=1}^3\left(\prod _{j=1}^k\frac{{\mu }_j+{\alpha }_j}{({\mu }_j+{\alpha }_j)+W_D^{\diamond }}\right)\frac{{\overline{S}}_k}{\sum _{j\in \mathcal{S}\setminus \{E\}}{\overline{S}}_j}{\mathcal{T}}_{0k}^{T\to D}{\overline{S}}_D\end{array}$(67) ), and multiplying (80(80) $\begin{array}{r}1=\frac{S_D^{\diamond }+ϵR_D^{\diamond }}{N_D^{\ast }}\frac{f(K{\mathcal{T}}_0)}{{\mathcal{T}}_0^{T\to D}}{\mathcal{T}}_0.\end{array}$(80) ) by K it follows that (81) $\begin{array}{r}K=\frac{S_D^{\diamond }+ϵR_D^{\diamond }}{N_D^{\ast }}\frac{P_T^{\ast }}{{\mathcal{T}}_0^{T\to D}}\mathrm{and}{W}_D^{\diamond }=\frac{S_D^{\diamond }+ϵR_D^{\diamond }}{N_D^{\ast }}\frac{P_T^{\ast }}{{\mathcal{T}}_0^{T\to D}}{\mathcal{T}}_0.\end{array}$(81) Therefore, our definition of $R_D^{\diamond }$ and ${\mathit{i}}_D^{\diamond }$ respectively in (74a(74a) $\begin{array}{rl}{\displaystyle R_D^{\diamond }({\mathcal{T}}_0,K)}& :={\mathcal{T}}_0^{T\to D}\frac{K{N}_D^{\ast }(K)}{ϵK{\mathcal{T}}_{0}f(K{\mathcal{T}}_0)+{\mu }_D{N}_D^{\ast }(K)}{\int }_0^{\mathrm{\infty }}\gamma (a){\pi }_D(a)\mathrm{d}a\end{array}$(74a) ) and (74c(74c) $\begin{array}{rl}{\displaystyle {\mathit{i}}_D^{\diamond }(a)}& :=K{\mathcal{T}}_0^{T\to D}{\pi }_D(a){\mathbf{v}}_1+K{\mathcal{T}}_0^{T\to D}{\int }_0^{\mathrm{\infty }}{\nu }_D(a){\pi }_D(a)\mathrm{d}a{\pi }_C(a){\mathbf{v}}_2.\end{array}$(74c) ) together with (81(81) $\begin{array}{r}K=\frac{S_D^{\diamond }+ϵR_D^{\diamond }}{N_D^{\ast }}\frac{P_T^{\ast }}{{\mathcal{T}}_0^{T\to D}}\mathrm{and}{W}_D^{\diamond }=\frac{S_D^{\diamond }+ϵR_D^{\diamond }}{N_D^{\ast }}\frac{P_T^{\ast }}{{\mathcal{T}}_0^{T\to D}}{\mathcal{T}}_0.\end{array}$(81) ) give us (82) $\begin{array}{r}\{\begin{array}{l}{\displaystyle R_D^{\diamond }=\frac{S_D^{\diamond }+ϵR_D^{\diamond }}{N_D^{\ast }}{P}_T^{\ast }\frac{N_D^{\ast }}{ϵP_T^{\ast }+{\mu }_D{N}_D^{\ast }}{\int }_0^{\mathrm{\infty }}\gamma (a){\pi }_D(a)\mathrm{d}a}\\ {\displaystyle {\mathit{i}}_D^{\diamond }(a)=\frac{S_D^{\diamond }+ϵR_D^{\diamond }}{N_D^{\ast }}{P}_T^{\ast }{\pi }_D(a){\mathbf{v}}_1+\frac{S_D^{\diamond }+ϵR_D^{\diamond }}{N_D^{\ast }}{P}_T^{\ast }{\int }_0^{\mathrm{\infty }}{\nu }_D(a){\pi }_D(a)\mathrm{d}a{\pi }_C(a){\mathbf{v}}_2.}\end{array}\end{array}$(82) Moreover, using (81(81) $\begin{array}{r}K=\frac{S_D^{\diamond }+ϵR_D^{\diamond }}{N_D^{\ast }}\frac{P_T^{\ast }}{{\mathcal{T}}_0^{T\to D}}\mathrm{and}{W}_D^{\diamond }=\frac{S_D^{\diamond }+ϵR_D^{\diamond }}{N_D^{\ast }}\frac{P_T^{\ast }}{{\mathcal{T}}_0^{T\to D}}{\mathcal{T}}_0.\end{array}$(81) ) one also note that our definition of $N_D^{\ast }$ in (73(73) $\begin{array}{r}{N}_D^{\ast }(K):=\frac{{\theta }_D}{{\mu }_D+K{\mathcal{T}}_0^{T\to D}({\int }_0^{\mathrm{\infty }}{\delta }_D(a){\pi }_D(a)\mathrm{d}a+{\int }_0^{\mathrm{\infty }}{\delta }_C(a){\pi }_C(a)\mathrm{d}a)}.\end{array}$(73) ) leads to (83) $\begin{array}{rl}{\mu }_D{N}_D^{\ast }& ={\theta }_D-(S_D^{\diamond }+ϵR_D^{\diamond })P_T^{\ast }{\int }_0^{\mathrm{\infty }}{\delta }_D(a){\pi }_D(a)\mathrm{d}a\\ & -(S_D^{\diamond }+ϵR_D^{\diamond })P_T^{\ast }{\int }_0^{\mathrm{\infty }}{\delta }_C(a){\pi }_C(a)\mathrm{d}a.\end{array}$(83) Next, we observe that the $R_D^{\diamond }$-equation in (82(82) $\begin{array}{r}\{\begin{array}{l}{\displaystyle R_D^{\diamond }=\frac{S_D^{\diamond }+ϵR_D^{\diamond }}{N_D^{\ast }}{P}_T^{\ast }\frac{N_D^{\ast }}{ϵP_T^{\ast }+{\mu }_D{N}_D^{\ast }}{\int }_0^{\mathrm{\infty }}\gamma (a){\pi }_D(a)\mathrm{d}a}\\ {\displaystyle {\mathit{i}}_D^{\diamond }(a)=\frac{S_D^{\diamond }+ϵR_D^{\diamond }}{N_D^{\ast }}{P}_T^{\ast }{\pi }_D(a){\mathbf{v}}_1+\frac{S_D^{\diamond }+ϵR_D^{\diamond }}{N_D^{\ast }}{P}_T^{\ast }{\int }_0^{\mathrm{\infty }}{\nu }_D(a){\pi }_D(a)\mathrm{d}a{\pi }_C(a){\mathbf{v}}_2.}\end{array}\end{array}$(82) ) is equivalent to (84) $\begin{array}{r}{\mu }_D{N}_D^{\ast }{R}_D^{\diamond }=-ϵP_T^{\ast }{R}_D^{\diamond }+(S_D^{\diamond }+ϵR_D^{\diamond })P_T^{\ast }{\int }_0^{\mathrm{\infty }}\gamma (a){\pi }_D(a)\mathrm{d}a.\end{array}$(84) Note that from the above formulas of $R_D^{\diamond }$ and ${\mathit{i}}_D^{\diamond }$ in (82(82) $\begin{array}{r}\{\begin{array}{l}{\displaystyle R_D^{\diamond }=\frac{S_D^{\diamond }+ϵR_D^{\diamond }}{N_D^{\ast }}{P}_T^{\ast }\frac{N_D^{\ast }}{ϵP_T^{\ast }+{\mu }_D{N}_D^{\ast }}{\int }_0^{\mathrm{\infty }}\gamma (a){\pi }_D(a)\mathrm{d}a}\\ {\displaystyle {\mathit{i}}_D^{\diamond }(a)=\frac{S_D^{\diamond }+ϵR_D^{\diamond }}{N_D^{\ast }}{P}_T^{\ast }{\pi }_D(a){\mathbf{v}}_1+\frac{S_D^{\diamond }+ϵR_D^{\diamond }}{N_D^{\ast }}{P}_T^{\ast }{\int }_0^{\mathrm{\infty }}{\nu }_D(a){\pi }_D(a)\mathrm{d}a{\pi }_C(a){\mathbf{v}}_2.}\end{array}\end{array}$(82) ) we have the following identities (85) $\begin{array}{r}\{\begin{array}{l}{\displaystyle {\int }_0^{\mathrm{\infty }}{\nu }_D(a)⟨{\mathbf{v}}_1,{\mathit{i}}_D^{\diamond }(a)⟩\mathrm{d}a=\frac{S_D^{\diamond }+ϵR_D^{\diamond }}{N_D^{\ast }}{P}_T^{\ast }{\int }_0^{\mathrm{\infty }}{\nu }_D(a){\pi }_D(a)\mathrm{d}a}\\ {\displaystyle {\int }_0^{\mathrm{\infty }}\gamma (a)⟨{\mathbf{v}}_1,{\mathit{i}}_D^{\diamond }(a)⟩\mathrm{d}a=\frac{S_D^{\diamond }+ϵR_D^{\diamond }}{N_D^{\ast }}{P}_T^{\ast }{\int }_0^{\mathrm{\infty }}\gamma (a){\pi }_D(a)\mathrm{d}a}\\ {\displaystyle {\int }_0^{\mathrm{\infty }}{\delta }_D(a)⟨{\mathbf{v}}_1,{\mathit{i}}_D^{\diamond }(a)⟩\mathrm{d}a=\frac{S_D^{\diamond }+ϵR_D^{\diamond }}{N_D^{\ast }}{P}_T^{\ast }{\int }_0^{\mathrm{\infty }}{\delta }_D(a){\pi }_D(a)\mathrm{d}a}\\ {\displaystyle {\int }_0^{\mathrm{\infty }}{\delta }_C(a)⟨{\mathbf{v}}_2,{\mathit{i}}_D^{\diamond }(a)⟩\mathrm{d}a=\frac{S_D^{\diamond }+ϵR_D^{\diamond }}{N_D^{\ast }}{P}_T^{\ast }{\int }_0^{\mathrm{\infty }}{\nu }_D(a){\pi }_D(a)\mathrm{d}a}\\ {\displaystyle \times {\int }_0^{\mathrm{\infty }}{\delta }_C(a){\pi }_C(a)\mathrm{d}a}\end{array}\end{array}$(85) from where we can rewrite (83(83) $\begin{array}{rl}{\mu }_D{N}_D^{\ast }& ={\theta }_D-(S_D^{\diamond }+ϵR_D^{\diamond })P_T^{\ast }{\int }_0^{\mathrm{\infty }}{\delta }_D(a){\pi }_D(a)\mathrm{d}a\\ & -(S_D^{\diamond }+ϵR_D^{\diamond })P_T^{\ast }{\int }_0^{\mathrm{\infty }}{\delta }_C(a){\pi }_C(a)\mathrm{d}a.\end{array}$(83) ), (84(84) $\begin{array}{r}{\mu }_D{N}_D^{\ast }{R}_D^{\diamond }=-ϵP_T^{\ast }{R}_D^{\diamond }+(S_D^{\diamond }+ϵR_D^{\diamond })P_T^{\ast }{\int }_0^{\mathrm{\infty }}\gamma (a){\pi }_D(a)\mathrm{d}a.\end{array}$(84) ) and the ${\mathit{i}}_D^{\ast }$-equation in (82(82) $\begin{array}{r}\{\begin{array}{l}{\displaystyle R_D^{\diamond }=\frac{S_D^{\diamond }+ϵR_D^{\diamond }}{N_D^{\ast }}{P}_T^{\ast }\frac{N_D^{\ast }}{ϵP_T^{\ast }+{\mu }_D{N}_D^{\ast }}{\int }_0^{\mathrm{\infty }}\gamma (a){\pi }_D(a)\mathrm{d}a}\\ {\displaystyle {\mathit{i}}_D^{\diamond }(a)=\frac{S_D^{\diamond }+ϵR_D^{\diamond }}{N_D^{\ast }}{P}_T^{\ast }{\pi }_D(a){\mathbf{v}}_1+\frac{S_D^{\diamond }+ϵR_D^{\diamond }}{N_D^{\ast }}{P}_T^{\ast }{\int }_0^{\mathrm{\infty }}{\nu }_D(a){\pi }_D(a)\mathrm{d}a{\pi }_C(a){\mathbf{v}}_2.}\end{array}\end{array}$(82) ) as follow (86) $\begin{array}{r}\{\begin{array}{l}{\displaystyle {\mu }_D{N}_D^{\ast }={\displaystyle {\theta }_D-N_D^{\ast }{\displaystyle {\int }_0^{\mathrm{\infty }}{\delta }_D(a)⟨{\mathbf{v}}_1,{\mathit{i}}_D^{\diamond }(a)⟩\mathrm{d}a-N_D^{\ast }{\displaystyle {\int }_0^{\mathrm{\infty }}{\delta }_C(a)⟨{\mathbf{v}}_2,{\mathit{i}}_D^{\diamond }(a)⟩\mathrm{d}a}}}}\\ {\mu }_D{N}_D^{\ast }{R}_D^{\diamond }={\displaystyle -ϵP_T^{\ast }{R}_D^{\diamond }+N_D^{\ast }{\int }_0^{\mathrm{\infty }}\gamma (a)⟨{\mathbf{v}}_1,{\mathit{i}}_D^{\diamond }(a)⟩\mathrm{d}a}\\ {\mathit{i}}_D^{\diamond }(a)={\displaystyle \frac{S_D^{\diamond }+ϵR_D^{\diamond }}{N_D^{\ast }}{P}_T^{\ast }{\pi }_D(a){\mathbf{v}}_1+{\int }_0^{\mathrm{\infty }}{\nu }_D(a)⟨{\mathbf{v}}_1,{\mathit{i}}_D^{\diamond }(a)⟩\mathrm{d}a{\pi }_C(a){\mathbf{v}}_2.}\end{array}\end{array}$(86) Next, observe that (87) $\begin{array}{r}\{\begin{array}{l}{\mathit{i}}_D^{\diamond }(a{)}^{\prime }=-({\mu }_D+\mathit{\theta }(a)){\mathit{i}}_D^{\diamond }(a)\\ {\displaystyle {\mathit{i}}_D^{\diamond }(0)=\frac{S_D^{\diamond }+ϵR_D^{\diamond }}{N_D^{\ast }}{P}_T^{\ast }{\mathbf{v}}_1+{\int }_0^{\mathrm{\infty }}{\nu }_D(a)⟨{\mathbf{v}}_1,{\mathit{i}}_D^{\diamond }(a)⟩\mathrm{d}a{\mathbf{v}}_2}\end{array}\end{array}$(87) so that by integrating (87(87) $\begin{array}{r}\{\begin{array}{l}{\mathit{i}}_D^{\diamond }(a{)}^{\prime }=-({\mu }_D+\mathit{\theta }(a)){\mathit{i}}_D^{\diamond }(a)\\ {\displaystyle {\mathit{i}}_D^{\diamond }(0)=\frac{S_D^{\diamond }+ϵR_D^{\diamond }}{N_D^{\ast }}{P}_T^{\ast }{\mathbf{v}}_1+{\int }_0^{\mathrm{\infty }}{\nu }_D(a)⟨{\mathbf{v}}_1,{\mathit{i}}_D^{\diamond }(a)⟩\mathrm{d}a{\mathbf{v}}_2}\end{array}\end{array}$(87) ) from 0 to $+\mathrm{\infty }$ we obtain $\begin{array}{r}-\frac{S_D^{\diamond }+ϵR_D^{\diamond }}{N_D^{\ast }}{P}_T^{\ast }{\mathbf{v}}_1-{\int }_0^{\mathrm{\infty }}{\nu }_D(a)⟨{\mathbf{v}}_1,{\mathit{i}}_D^{\diamond }(a)⟩\mathrm{d}a{\mathbf{v}}_2=-{\int }_0^{\mathrm{\infty }}({\mu }_D+\mathit{\theta }(a)){\mathit{i}}_D^{\diamond }(a)\mathrm{d}a\end{array}$and by using (85(85) $\begin{array}{r}\{\begin{array}{l}{\displaystyle {\int }_0^{\mathrm{\infty }}{\nu }_D(a)⟨{\mathbf{v}}_1,{\mathit{i}}_D^{\diamond }(a)⟩\mathrm{d}a=\frac{S_D^{\diamond }+ϵR_D^{\diamond }}{N_D^{\ast }}{P}_T^{\ast }{\int }_0^{\mathrm{\infty }}{\nu }_D(a){\pi }_D(a)\mathrm{d}a}\\ {\displaystyle {\int }_0^{\mathrm{\infty }}\gamma (a)⟨{\mathbf{v}}_1,{\mathit{i}}_D^{\diamond }(a)⟩\mathrm{d}a=\frac{S_D^{\diamond }+ϵR_D^{\diamond }}{N_D^{\ast }}{P}_T^{\ast }{\int }_0^{\mathrm{\infty }}\gamma (a){\pi }_D(a)\mathrm{d}a}\\ {\displaystyle {\int }_0^{\mathrm{\infty }}{\delta }_D(a)⟨{\mathbf{v}}_1,{\mathit{i}}_D^{\diamond }(a)⟩\mathrm{d}a=\frac{S_D^{\diamond }+ϵR_D^{\diamond }}{N_D^{\ast }}{P}_T^{\ast }{\int }_0^{\mathrm{\infty }}{\delta }_D(a){\pi }_D(a)\mathrm{d}a}\\ {\displaystyle {\int }_0^{\mathrm{\infty }}{\delta }_C(a)⟨{\mathbf{v}}_2,{\mathit{i}}_D^{\diamond }(a)⟩\mathrm{d}a=\frac{S_D^{\diamond }+ϵR_D^{\diamond }}{N_D^{\ast }}{P}_T^{\ast }{\int }_0^{\mathrm{\infty }}{\nu }_D(a){\pi }_D(a)\mathrm{d}a}\\ {\displaystyle \times {\int }_0^{\mathrm{\infty }}{\delta }_C(a){\pi }_C(a)\mathrm{d}a}\end{array}\end{array}$(85) ) it follows that (88) $\begin{array}{rl}{\mu }_D{\int }_0^{\mathrm{\infty }}⟨\mathbf{1},{\mathit{i}}_D^{\diamond }(a)⟩\mathrm{d}a& ={\displaystyle \frac{S_D^{\diamond }+ϵR_D^{\diamond }}{N_D^{\ast }}{P}_T^{\ast }-{\int }_0^{\mathrm{\infty }}(\gamma (a)+{\delta }_D(a))⟨{\mathbf{v}}_1,{\mathit{i}}_D^{\diamond }(a)⟩\mathrm{d}a}\\ & -{\int }_0^{\mathrm{\infty }}{\delta }_C(a)⟨{\mathbf{v}}_2,{\mathit{i}}_D^{\diamond }(a)⟩\mathrm{d}a.\end{array}$(88) Next, multiply (89(89) $\begin{array}{r}{\mu }_D{N}_D^{\ast }{\int }_0^{\mathrm{\infty }}⟨\mathbf{1},{\mathit{i}}_D^{\diamond }(a)⟩\mathrm{d}a=(S_D^{\diamond }+ϵR_D^{\diamond })P_T^{\ast }-{\theta }_D+{\mu }_D{N}_D^{\ast }-N_D^{\ast }{\int }_0^{\mathrm{\infty }}\gamma (a)⟨{\mathbf{v}}_1,{\mathit{i}}_D^{\diamond }(a)⟩\mathrm{d}a\end{array}$(89) ) by $N_D^{\ast }$ and use the first equation of (86(86) $\begin{array}{r}\{\begin{array}{l}{\displaystyle {\mu }_D{N}_D^{\ast }={\displaystyle {\theta }_D-N_D^{\ast }{\displaystyle {\int }_0^{\mathrm{\infty }}{\delta }_D(a)⟨{\mathbf{v}}_1,{\mathit{i}}_D^{\diamond }(a)⟩\mathrm{d}a-N_D^{\ast }{\displaystyle {\int }_0^{\mathrm{\infty }}{\delta }_C(a)⟨{\mathbf{v}}_2,{\mathit{i}}_D^{\diamond }(a)⟩\mathrm{d}a}}}}\\ {\mu }_D{N}_D^{\ast }{R}_D^{\diamond }={\displaystyle -ϵP_T^{\ast }{R}_D^{\diamond }+N_D^{\ast }{\int }_0^{\mathrm{\infty }}\gamma (a)⟨{\mathbf{v}}_1,{\mathit{i}}_D^{\diamond }(a)⟩\mathrm{d}a}\\ {\mathit{i}}_D^{\diamond }(a)={\displaystyle \frac{S_D^{\diamond }+ϵR_D^{\diamond }}{N_D^{\ast }}{P}_T^{\ast }{\pi }_D(a){\mathbf{v}}_1+{\int }_0^{\mathrm{\infty }}{\nu }_D(a)⟨{\mathbf{v}}_1,{\mathit{i}}_D^{\diamond }(a)⟩\mathrm{d}a{\pi }_C(a){\mathbf{v}}_2.}\end{array}\end{array}$(86) ) to obtain (89) $\begin{array}{r}{\mu }_D{N}_D^{\ast }{\int }_0^{\mathrm{\infty }}⟨\mathbf{1},{\mathit{i}}_D^{\diamond }(a)⟩\mathrm{d}a=(S_D^{\diamond }+ϵR_D^{\diamond })P_T^{\ast }-{\theta }_D+{\mu }_D{N}_D^{\ast }-N_D^{\ast }{\int }_0^{\mathrm{\infty }}\gamma (a)⟨{\mathbf{v}}_1,{\mathit{i}}_D^{\diamond }(a)⟩\mathrm{d}a\end{array}$(89) Hence, summing (89(89) $\begin{array}{r}{\mu }_D{N}_D^{\ast }{\int }_0^{\mathrm{\infty }}⟨\mathbf{1},{\mathit{i}}_D^{\diamond }(a)⟩\mathrm{d}a=(S_D^{\diamond }+ϵR_D^{\diamond })P_T^{\ast }-{\theta }_D+{\mu }_D{N}_D^{\ast }-N_D^{\ast }{\int }_0^{\mathrm{\infty }}\gamma (a)⟨{\mathbf{v}}_1,{\mathit{i}}_D^{\diamond }(a)⟩\mathrm{d}a\end{array}$(89) ) and the second equation of (86(86) $\begin{array}{r}\{\begin{array}{l}{\displaystyle {\mu }_D{N}_D^{\ast }={\displaystyle {\theta }_D-N_D^{\ast }{\displaystyle {\int }_0^{\mathrm{\infty }}{\delta }_D(a)⟨{\mathbf{v}}_1,{\mathit{i}}_D^{\diamond }(a)⟩\mathrm{d}a-N_D^{\ast }{\displaystyle {\int }_0^{\mathrm{\infty }}{\delta }_C(a)⟨{\mathbf{v}}_2,{\mathit{i}}_D^{\diamond }(a)⟩\mathrm{d}a}}}}\\ {\mu }_D{N}_D^{\ast }{R}_D^{\diamond }={\displaystyle -ϵP_T^{\ast }{R}_D^{\diamond }+N_D^{\ast }{\int }_0^{\mathrm{\infty }}\gamma (a)⟨{\mathbf{v}}_1,{\mathit{i}}_D^{\diamond }(a)⟩\mathrm{d}a}\\ {\mathit{i}}_D^{\diamond }(a)={\displaystyle \frac{S_D^{\diamond }+ϵR_D^{\diamond }}{N_D^{\ast }}{P}_T^{\ast }{\pi }_D(a){\mathbf{v}}_1+{\int }_0^{\mathrm{\infty }}{\nu }_D(a)⟨{\mathbf{v}}_1,{\mathit{i}}_D^{\diamond }(a)⟩\mathrm{d}a{\pi }_C(a){\mathbf{v}}_2.}\end{array}\end{array}$(86) ) it follows that (90) $\begin{array}{rl}{S}_D^{\diamond }{P}_T^{\ast }-{\theta }_D+{\mu }_D{N}_D^{\ast }& ={\displaystyle {\mu }_D{N}_D^{\ast }{\int }_0^{\mathrm{\infty }}⟨\mathbf{1},{\mathit{i}}_D^{\diamond }(a)⟩\mathrm{d}a+{\mu }_D{N}_D^{\ast }{R}_D^{\diamond }}\\ & ={\mu }_D{N}_D^{\ast }(1-S_D^{\diamond })\end{array}$(90) that is (91) $\begin{array}{r}0={\theta }_D-S_D^{\diamond }{P}_T^{\ast }-{\mu }_D{N}_D^{\ast }{S}_D^{\diamond }⟺S_D^{\diamond }=\frac{{\theta }_D}{{\mu }_D{N}_D^{\ast }+P_T^{\ast }}.\end{array}$(91) Finally, we infer from the ${\mathit{i}}_D^{\diamond }$-equation of (86(86) $\begin{array}{r}\{\begin{array}{l}{\displaystyle {\mu }_D{N}_D^{\ast }={\displaystyle {\theta }_D-N_D^{\ast }{\displaystyle {\int }_0^{\mathrm{\infty }}{\delta }_D(a)⟨{\mathbf{v}}_1,{\mathit{i}}_D^{\diamond }(a)⟩\mathrm{d}a-N_D^{\ast }{\displaystyle {\int }_0^{\mathrm{\infty }}{\delta }_C(a)⟨{\mathbf{v}}_2,{\mathit{i}}_D^{\diamond }(a)⟩\mathrm{d}a}}}}\\ {\mu }_D{N}_D^{\ast }{R}_D^{\diamond }={\displaystyle -ϵP_T^{\ast }{R}_D^{\diamond }+N_D^{\ast }{\int }_0^{\mathrm{\infty }}\gamma (a)⟨{\mathbf{v}}_1,{\mathit{i}}_D^{\diamond }(a)⟩\mathrm{d}a}\\ {\mathit{i}}_D^{\diamond }(a)={\displaystyle \frac{S_D^{\diamond }+ϵR_D^{\diamond }}{N_D^{\ast }}{P}_T^{\ast }{\pi }_D(a){\mathbf{v}}_1+{\int }_0^{\mathrm{\infty }}{\nu }_D(a)⟨{\mathbf{v}}_1,{\mathit{i}}_D^{\diamond }(a)⟩\mathrm{d}a{\pi }_C(a){\mathbf{v}}_2.}\end{array}\end{array}$(86) ), the equality (91(91) $\begin{array}{r}0={\theta }_D-S_D^{\diamond }{P}_T^{\ast }-{\mu }_D{N}_D^{\ast }{S}_D^{\diamond }⟺S_D^{\diamond }=\frac{{\theta }_D}{{\mu }_D{N}_D^{\ast }+P_T^{\ast }}.\end{array}$(91) ) together with our definition of $R_D^{\diamond }$, and $N_D^{\ast }$ that (62(62) $\begin{array}{r}\{\begin{array}{l}{\displaystyle N_D^{\ast }=\frac{{\theta }_D}{{\mu }_D+U_D^{\diamond }+U_C^{\diamond }}}\\ {\displaystyle S_D^{\diamond }=\frac{{\theta }_D}{{\mu }_D{N}_D^{\ast }+P_T^{\ast }}}\\ {\displaystyle R_D^{\diamond }=\frac{N_D^{\ast }{P}_D^{\diamond }}{ϵP_T^{\ast }+{\mu }_D{N}_D^{\ast }}}\\ {\displaystyle {\mathit{i}}_D^{\diamond }(a)=\frac{S_D^{\diamond }+ϵR_D^{\diamond }}{N_D^{\ast }}{P}_T^{\ast }{\pi }_D(a){\mathbf{v}}_1+Q_D^{\diamond }{\pi }_C(a){\mathbf{v}}_2.}\end{array}\end{array}$(62) ) is satisfied. The proof is completed by using the fact that we have set $P_T^{\ast }=W_D^{\diamond }f(W_D^{\diamond })$ with $W_D^{\diamond }=K{\mathcal{T}}_0$ and f given by (67(67) $\begin{array}{r}f(W_D^{\diamond }):=\sum _{k=1}^3\left(\prod _{j=1}^k\frac{{\mu }_j+{\alpha }_j}{({\mu }_j+{\alpha }_j)+W_D^{\diamond }}\right)\frac{{\overline{S}}_k}{\sum _{j\in \mathcal{S}\setminus \{E\}}{\overline{S}}_j}{\mathcal{T}}_{0k}^{T\to D}{\overline{S}}_D\end{array}$(67) ).

Thanks to Lemma 8.2, the existence of positive equilibrium to (7(7) $\begin{array}{r}\{\begin{array}{l}{\displaystyle \frac{\mathrm{d}{S}_D(t)}{\mathrm{d}t}={\theta }_D-\frac{S_D(t)}{N_D(t)}{\int }_0^{\mathrm{\infty }}{\beta }_T(a)⟨\mathbf{1},{\mathit{i}}_T(t,a)⟩\mathrm{d}a-{\mu }_D{S}_D(t),}\\ {\displaystyle \frac{\mathrm{\partial }{\mathit{i}}_D(t,a)}{\mathrm{\partial }t}+\frac{\mathrm{\partial }{\mathit{i}}_D(t,a)}{\mathrm{\partial }a}=-({\mu }_D+\mathit{\theta }(a)){\mathit{i}}_D(t,a),}\\ {\mathit{i}}_D(t,0)={\displaystyle \frac{S_D(t)+ϵR_D(t)}{N_D(t)}{\int }_0^{\mathrm{\infty }}{\beta }_T(a)⟨\mathbf{1},{\mathit{i}}_T(t,a)⟩\mathrm{d}a{\mathbf{v}}_1}\\ {\displaystyle +{\int }_0^{\mathrm{\infty }}{\nu }_D(a)⟨{\mathbf{v}}_1,{\mathit{i}}_D(t,a)⟩\mathrm{d}a{\mathbf{v}}_2,}\\ {\displaystyle \frac{\mathrm{d}{R}_D(t)}{\mathrm{d}t}={\int }_0^{\mathrm{\infty }}\gamma (a)⟨{\mathbf{v}}_1,{\mathit{i}}_D(t,a)⟩\mathrm{d}a-ϵ\frac{R_D(t)}{N_D(t)}{\int }_0^{\mathrm{\infty }}{\beta }_T(a)⟨\mathbf{1},{\mathit{i}}_T(t,a)⟩\mathrm{d}a-{\mu }_D{R}_D(t),}\end{array}\end{array}$(7) )–(8(8) $\begin{array}{r}\{\begin{array}{l}{\displaystyle \frac{\mathrm{d}{S}_{TE}(t)}{\mathrm{d}t}=r_E(1-\frac{S_{TE}(t)}{K}){\displaystyle (⟨{\mathit{S}}_T(t),{\mathbf{e}}_3⟩+{\int }_0^{\mathrm{\infty }}⟨{\mathit{i}}_T(t,a),{\mathbf{e}}_3⟩\mathrm{d}a)-({\alpha }_E+{\mu }_E)S_{TE}(t),}}\\ {\displaystyle \frac{\mathrm{d}{\mathit{S}}_T(t)}{\mathrm{d}t}={\alpha }_E{S}_{TE}(t){\mathbf{e}}_1-\mathbf{M}{\mathit{S}}_T(t)-\frac{{\mathit{S}}_T(t)}{N_D(t)}{\int }_0^{\mathrm{\infty }}⟨{\mathit{\beta }}_D(a),{\mathit{i}}_D(t,a)⟩\mathrm{d}a,}\\ {\displaystyle \frac{\mathrm{\partial }{\mathit{i}}_T(t,a)}{\mathrm{\partial }t}+\frac{\mathrm{\partial }{\mathit{i}}_T(t,a)}{\mathrm{\partial }a}=-\mathbf{M}{\mathit{i}}_T(t,a),}\\ {\mathit{i}}_T(t,0)={\displaystyle \frac{{\mathit{S}}_T(t)}{N_D(t)}{\int }_0^{\mathrm{\infty }}⟨{\mathit{\beta }}_D(a),{\mathit{i}}_D(t,a)⟩\mathrm{d}a,}\end{array}\end{array}$(8) ) is subjected to the existence of K>0 such that (92) $\begin{array}{r}\mathrm{\Delta }({\mathcal{T}}_0,K)=1\end{array}$(92) where $\mathrm{\Delta }({\mathcal{T}}_0,K)$ is given by (76(76) $\begin{array}{r}\mathrm{\Delta }({\mathcal{T}}_0,K):=\frac{1-(1-ϵ)R_D^{\diamond }({\mathcal{T}}_0,K)-I_D^{\diamond }({\mathcal{T}}_0,K)}{N_D^{\ast }(K)}\frac{f(K{\mathcal{T}}_0)}{{\mathcal{T}}_0^{T\to D}}{\mathcal{T}}_0\end{array}$(76) ) and (77(77) $\begin{array}{r}{I}_D^{\diamond }({\mathcal{T}}_0,K):=K{\mathcal{T}}_0^{T\to D}{\int }_0^{\mathrm{\infty }}{\pi }_D(a)\mathrm{d}a+K{\mathcal{T}}_0^{T\to D}{\int }_0^{\mathrm{\infty }}{\nu }_D(a){\pi }_D(a)\mathrm{d}a{\int }_0^{\mathrm{\infty }}{\pi }_C(a)\mathrm{d}a.\end{array}$(77) ). Before studying the existence of solution to (92(92) $\begin{array}{r}\mathrm{\Delta }({\mathcal{T}}_0,K)=1\end{array}$(92) ), we first observe that from (73(73) $\begin{array}{r}{N}_D^{\ast }(K):=\frac{{\theta }_D}{{\mu }_D+K{\mathcal{T}}_0^{T\to D}({\int }_0^{\mathrm{\infty }}{\delta }_D(a){\pi }_D(a)\mathrm{d}a+{\int }_0^{\mathrm{\infty }}{\delta }_C(a){\pi }_C(a)\mathrm{d}a)}.\end{array}$(73) ), (74a(74a) $\begin{array}{rl}{\displaystyle R_D^{\diamond }({\mathcal{T}}_0,K)}& :={\mathcal{T}}_0^{T\to D}\frac{K{N}_D^{\ast }(K)}{ϵK{\mathcal{T}}_{0}f(K{\mathcal{T}}_0)+{\mu }_D{N}_D^{\ast }(K)}{\int }_0^{\mathrm{\infty }}\gamma (a){\pi }_D(a)\mathrm{d}a\end{array}$(74a) ) and (77(77) $\begin{array}{r}{I}_D^{\diamond }({\mathcal{T}}_0,K):=K{\mathcal{T}}_0^{T\to D}{\int }_0^{\mathrm{\infty }}{\pi }_D(a)\mathrm{d}a+K{\mathcal{T}}_0^{T\to D}{\int }_0^{\mathrm{\infty }}{\nu }_D(a){\pi }_D(a)\mathrm{d}a{\int }_0^{\mathrm{\infty }}{\pi }_C(a)\mathrm{d}a.\end{array}$(77) ) we have (93) $\begin{array}{r}\{\begin{array}{l}{N}_D^{\ast }(0)=\frac{{\theta }_D}{{\mu }_D}={\overline{S}}_D\\ R_D^{\diamond }({\mathcal{T}}_0,0)=I_D^{\diamond }({\mathcal{T}}_0,0)=0,\mathrm{\forall }{\mathcal{T}}_0\ge 0.\end{array}\end{array}$(93) We also note that for each ${\mathcal{T}}_0>0$ we have (94) $\begin{array}{r}\underset{K\to +\mathrm{\infty }}{lim}{R}_D^{\diamond }({\mathcal{T}}_0,K)=l>0\mathrm{and}\underset{K\to +\mathrm{\infty }}{lim}{I}_D^{\diamond }({\mathcal{T}}_0,K)=+\mathrm{\infty }\end{array}$(94) with $\begin{array}{r}l:=\frac{1}{ϵ{\mathcal{T}}_0^{T\to D}}\frac{{\theta }_D{\int }_0^{\mathrm{\infty }}\gamma (a){\pi }_D(a)\mathrm{d}a}{{\int }_0^{\mathrm{\infty }}{\delta }_D(a){\pi }_D(a)\mathrm{d}a+{\int }_0^{\mathrm{\infty }}{\delta }_C(a){\pi }_C(a)\mathrm{d}a}.\end{array}$Therefore, using the equality ${\mathcal{T}}_0^{T\to D}=\frac{f(0)}{{\overline{S}}_D}$ it follows from (93(93) $\begin{array}{r}\{\begin{array}{l}{N}_D^{\ast }(0)=\frac{{\theta }_D}{{\mu }_D}={\overline{S}}_D\\ R_D^{\diamond }({\mathcal{T}}_0,0)=I_D^{\diamond }({\mathcal{T}}_0,0)=0,\mathrm{\forall }{\mathcal{T}}_0\ge 0.\end{array}\end{array}$(93) ) and (76(76) $\begin{array}{r}\mathrm{\Delta }({\mathcal{T}}_0,K):=\frac{1-(1-ϵ)R_D^{\diamond }({\mathcal{T}}_0,K)-I_D^{\diamond }({\mathcal{T}}_0,K)}{N_D^{\ast }(K)}\frac{f(K{\mathcal{T}}_0)}{{\mathcal{T}}_0^{T\to D}}{\mathcal{T}}_0\end{array}$(76) ) that (95) $\begin{array}{r}\mathrm{\Delta }({\mathcal{T}}_0,0)={\mathcal{T}}_0\mathrm{and}\underset{K\to +\mathrm{\infty }}{lim}\mathrm{\Delta }({\mathcal{T}}_0,K)=-\mathrm{\infty }.\end{array}$(95) It is now clear from (95(95) $\begin{array}{r}\mathrm{\Delta }({\mathcal{T}}_0,0)={\mathcal{T}}_0\mathrm{and}\underset{K\to +\mathrm{\infty }}{lim}\mathrm{\Delta }({\mathcal{T}}_0,K)=-\mathrm{\infty }.\end{array}$(95) ) that for each ${\mathcal{T}}_0>0$ there exists $\overline{K}:=\overline{K}({\mathcal{T}}_0)>0$ such that $\mathrm{\Delta }({\mathcal{T}}_0,\overline{K})=0$. In particular, if ${\mathcal{T}}_0>1$ then there exists $K_0\in (0,\overline{K})$ such that $\mathrm{\Delta }({\mathcal{T}}_0,K_0)=1$ leading to the existence of an endemic equilibrium.

On the forward and backward bifurcations: In the following, we deal with the existence of forward and backward bifurcation at ${\mathcal{T}}_0=1$. Roughly speaking, we have to prove that there exists a positive map $K_0$ defined in some right neighbourhood (forward bifurcation) or left neighbourhood (backward bifurcation) of ${\mathcal{T}}_0=1$ such that $\mathrm{\Delta }({\mathcal{T}}_0,K_0({\mathcal{T}}_0))=1$ and $K_0(1)=0$. Since $\mathrm{\Delta }(1,0)=1$, one can prove the existence of forward and backward bifurcations using implicit function theorem which is reduced here to the study of ${\mathrm{\partial }}_K\mathrm{\Delta }(1,0)$. This motivates the following lemma.

Lemma 8.3

The following property is satisfied (96) $\begin{array}{r}{\mathrm{\partial }}_{{\mathcal{T}}_0}\mathrm{\Delta }({\mathcal{T}}_0,0)=1,\mathrm{\forall }{\mathcal{T}}_0>0\end{array}$(96) and $\begin{array}{rl}{\mathrm{\partial }}_K\mathrm{\Delta }(1,0)& =\frac{1}{{\mu }_D}{\mathcal{T}}_0^{T\to D}({\int }_0^{\mathrm{\infty }}{\delta }_D(a){\pi }_D(a)\mathrm{d}a+{\int }_0^{\mathrm{\infty }}{\delta }_C(a){\pi }_C(a)\mathrm{d}a)\\ & -\sum _{k=1}^3\left(\sum _{j=1}^k\frac{1}{{\mu }_j+{\alpha }_j}\right)\frac{\frac{{\overline{S}}_k}{\sum _{j\in \mathcal{S}\setminus \{E\}}{\overline{S}}_j}{\mathcal{T}}_{0k}^{T\to D}}{{\mathcal{T}}_0^{T\to D}}\\ & -(1-ϵ)\frac{1}{{\mu }_D}{\mathcal{T}}_0^{T\to D}{\int }_0^{\mathrm{\infty }}\gamma (a){\pi }_D(a)\mathrm{d}a\\ & -{\mathcal{T}}_0^{T\to D}{\int }_0^{\mathrm{\infty }}{\pi }_D(a)\mathrm{d}a-{\mathcal{T}}_0^{T\to D}{\int }_0^{\mathrm{\infty }}{\nu }_D(a){\pi }_D(a)\mathrm{d}a{\int }_0^{\mathrm{\infty }}{\pi }_C(a)\mathrm{d}a.\end{array}$

Proof.

Since $\mathrm{\Delta }({\mathcal{T}}_0,0)={\mathcal{T}}_0$ for each ${\mathcal{T}}_0>0$ the first equality of the lemma follows. To compute ${\mathrm{\partial }}_K\mathrm{\Delta }(1,0)$, it is enough to take the derivative of $\mathrm{\Delta }(1,K)$ with respect to K at K = 0. To do so let us first note that by using (76(76) $\begin{array}{r}\mathrm{\Delta }({\mathcal{T}}_0,K):=\frac{1-(1-ϵ)R_D^{\diamond }({\mathcal{T}}_0,K)-I_D^{\diamond }({\mathcal{T}}_0,K)}{N_D^{\ast }(K)}\frac{f(K{\mathcal{T}}_0)}{{\mathcal{T}}_0^{T\to D}}{\mathcal{T}}_0\end{array}$(76) ) we have $\begin{array}{r}{N}_D^{\ast }(K)\mathrm{\Delta }(1,K)=\frac{f(K)}{{\mathcal{T}}_0^{T\to D}}-[(1-ϵ)R_D^{\diamond }(1,K)+I_D^{\diamond }(1,K)]\frac{f(K)}{{\mathcal{T}}_0^{T\to D}}.\end{array}$Recalling that $\mathrm{\Delta }(1,0)=1$, $N_D^{\ast }(0)=\frac{{\theta }_D}{{\mu }_D}={\overline{S}}_D$, $R_D^{\diamond }(1,0)=I_D^{\diamond }(1,0)=0$ and taking the derivative of $N_D^{\ast }(K)\mathrm{\Delta }(1,K)$ with respect to K at K = 0 it follows that $\begin{array}{r}{\overline{S}}_D{\mathrm{\partial }}_K\mathrm{\Delta }(1,0)=-\frac{\mathrm{d}{N}_D^{\ast }(0)}{\mathrm{d}K}+\frac{f^{\prime }(0)}{{\mathcal{T}}_0^{T\to D}}-{\overline{S}}_D[(1-ϵ){\mathrm{\partial }}_K{R}_D^{\diamond }(1,0)+{\mathrm{\partial }}_K{I}_D^{\diamond }(1,0)].\end{array}$By straightforward computations, we obtain from (73(73) $\begin{array}{r}{N}_D^{\ast }(K):=\frac{{\theta }_D}{{\mu }_D+K{\mathcal{T}}_0^{T\to D}({\int }_0^{\mathrm{\infty }}{\delta }_D(a){\pi }_D(a)\mathrm{d}a+{\int }_0^{\mathrm{\infty }}{\delta }_C(a){\pi }_C(a)\mathrm{d}a)}.\end{array}$(73) ) and (77(77) $\begin{array}{r}{I}_D^{\diamond }({\mathcal{T}}_0,K):=K{\mathcal{T}}_0^{T\to D}{\int }_0^{\mathrm{\infty }}{\pi }_D(a)\mathrm{d}a+K{\mathcal{T}}_0^{T\to D}{\int }_0^{\mathrm{\infty }}{\nu }_D(a){\pi }_D(a)\mathrm{d}a{\int }_0^{\mathrm{\infty }}{\pi }_C(a)\mathrm{d}a.\end{array}$(77) ) that $\begin{array}{r}\frac{\mathrm{d}{N}_D^{\ast }(0)}{\mathrm{d}K}=-\frac{{\overline{S}}_D}{{\mu }_D}{\mathcal{T}}_0^{T\to D}({\int }_0^{\mathrm{\infty }}{\delta }_D(a){\pi }_D(a)\mathrm{d}a+{\int }_0^{\mathrm{\infty }}{\delta }_C(a){\pi }_C(a)\mathrm{d}a)\end{array}$and $\begin{array}{r}{\mathrm{\partial }}_K{I}_D^{\diamond }(1,0)={\displaystyle {\mathcal{T}}_0^{T\to D}{\int }_0^{\mathrm{\infty }}{\pi }_D(a)\mathrm{d}a+{\mathcal{T}}_0^{T\to D}{\int }_0^{\mathrm{\infty }}{\nu }_D(a){\pi }_C(a)\mathrm{d}a{\int }_0^{\mathrm{\infty }}{\pi }_C(a)\mathrm{d}a.}\end{array}$To complete the proof it remains to compute ${\mathrm{\partial }}_K{R}_D^{\diamond }(1,0)$ and $f^{\prime }(0)$. We first compute ${\mathrm{\partial }}_K{R}_D^{\diamond }(1,0)$. To this end, we compute the derivative of $R_D^{\diamond }(1,K)$ from (74a(74a) $\begin{array}{rl}{\displaystyle R_D^{\diamond }({\mathcal{T}}_0,K)}& :={\mathcal{T}}_0^{T\to D}\frac{K{N}_D^{\ast }(K)}{ϵK{\mathcal{T}}_{0}f(K{\mathcal{T}}_0)+{\mu }_D{N}_D^{\ast }(K)}{\int }_0^{\mathrm{\infty }}\gamma (a){\pi }_D(a)\mathrm{d}a\end{array}$(74a) ) with respect to K and set K = 0 to obtain $\begin{array}{r}{\mathrm{\partial }}_K{R}_D^{\diamond }(1,0)=\frac{1}{{\mu }_D}{\mathcal{T}}_0^{T\to D}{\int }_0^{\mathrm{\infty }}\gamma (a){\pi }_D(a)\mathrm{d}a.\end{array}$To compute $f^{\prime }(0)$ recall that $\begin{array}{r}f(W_D^{\diamond }):=\sum _{k=1}^3\left(\prod _{j=1}^k\frac{{\mu }_j+{\alpha }_j}{({\mu }_j+{\alpha }_j)+W_D^{\diamond }}\right)\frac{{\overline{S}}_k}{\sum _{j\in \mathcal{S}\setminus \{E\}}{\overline{S}}_j}{\mathcal{T}}_{0k}^{T\to D}{\overline{S}}_D\end{array}$so that setting $\begin{array}{r}{\hat{f}}_k(W_D^{\diamond }):=\prod _{j=1}^k{g}_j(W_D^{\diamond }),k=1,2,3,g_j(W_D^{\diamond }):=\frac{{\mu }_j+{\alpha }_j}{({\mu }_j+{\alpha }_j)+W_D^{\diamond }}\end{array}$we obtain $\begin{array}{r}f(W_D^{\diamond }):=\sum _{k=1}^3{\hat{f}}_k(W_D^{\diamond })\frac{{\overline{S}}_k}{\sum _{j\in \mathcal{S}\setminus \{E\}}{\overline{S}}_j}{\mathcal{T}}_{0k}^{T\to D}{\overline{S}}_D.\end{array}$Therefore, computing the derivative of $\ln ({\hat{f}}_k(W_D^{\diamond }))$ one obtains $\begin{array}{r}{\hat{f}}_k^{\prime }(W_D^{\diamond })={\hat{f}}_k(W_D^{\diamond })\sum _{j=1}^k\frac{g_j^{\prime }(W_D^{\diamond })}{g_j(W_D^{\diamond })}\end{array}$and since ${\hat{f}}_k(0)=g_k(0)=1$ it follows that $\begin{array}{r}{\hat{f}}_k^{\prime }(0)=\sum _{j=1}^k{g}_j^{\prime }(0)=-\sum _{j=1}^k\frac{1}{{\mu }_j+{\alpha }_j}\end{array}$that is (97) $\begin{array}{r}{f}^{\prime }(0)=-\sum _{k=1}^3\left(\sum _{j=1}^k\frac{1}{{\mu }_j+{\alpha }_j}\right)\frac{{\overline{S}}_k}{\sum _{j\in \mathcal{S}\setminus \{E\}}{\overline{S}}_j}{\mathcal{T}}_{0k}^{T\to D}{\overline{S}}_D.\end{array}$(97) The proof is completed.

Thanks to Lemma 8.3 one knows that if ${\mathrm{\partial }}_K\mathrm{\Delta }(1,0)=C_{\mathrm{bif}}^{-1}\ne 0$ then by the implicit function theorem there exists $\xi \in (0,1)$ and a smooth map $K_0:(1-\xi ,1+\xi )\to \mathbb{R}$ such that (98) $\begin{array}{r}\mathrm{\Delta }({\mathcal{T}}_0,K_0({\mathcal{T}}_0))=1,\mathrm{\forall }{\mathcal{T}}_0\in (1-\xi ,1+\xi )\mathrm{and}{K}_0(1)=0.\end{array}$(98) Thus, up to reduce ξ, it follows that the sign of $K_0$ on $(1-\xi ,1+\xi )$ is given by the sign of $\frac{\mathrm{d}{K}_0({\mathcal{T}}_0)}{\mathrm{d}{\mathcal{T}}_0}$ at ${\mathcal{T}}_0=1$. Differentiating (98(98) $\begin{array}{r}\mathrm{\Delta }({\mathcal{T}}_0,K_0({\mathcal{T}}_0))=1,\mathrm{\forall }{\mathcal{T}}_0\in (1-\xi ,1+\xi )\mathrm{and}{K}_0(1)=0.\end{array}$(98) ) with respect to ${\mathcal{T}}_0$ and using Lemma 8.3 we obtain $\begin{array}{r}\frac{\mathrm{d}{K}_0(1)}{\mathrm{d}{\mathcal{T}}_0}=-\frac{{\mathrm{\partial }}_{{\mathcal{T}}_0}\mathrm{\Delta }(1,0)}{{\mathrm{\partial }}_K\mathrm{\Delta }(1,0)}=-C_{\mathrm{bif}}\end{array}$

so that

• If $C_{\mathrm{bif}}>0$ then, up to reduce ξ, we have $K_0({\mathcal{T}}_0)>0$ for ${\mathcal{T}}_0\in (1-\xi ,1)$ which corresponds to a backward bifurcation ;

• If $C_{\mathrm{bif}}<0$ then, up to reduce ξ, $K_0({\mathcal{T}}_0)>0$ for ${\mathcal{T}}_0\in (1,1+\xi )$ which corresponds to a forward bifurcation.

Disclosure statement

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

Additional information

Funding

FBA is supported by the National Science Foundation under the EPSCOR Track 2 grant number 1920946.