Derdei Bichara

Abderrahman Iggidr

Laura Smith

Multi-stage vector-borne zoonoses models: A global analysis

Keywords: Mathematics Subject Classification: 92D30, 34D23, 34D20, 34D40, 34A34 Vector-borne zoonoses, Stage progression, Multi-host, One Health, Amplification effect, Global stability, Nonlinear dynamical systems

A class of models that describes the interactions between multiple host species and an arthropod vector is formulated and its dynamics investigated. A host-vector disease model where the host's infection is structured into n stages is formulated and a complete global dynamics analysis is provided. The basic reproduction number acts as a sharp threshold, that is, the disease-free equilibrium is globally asymptotically stable (GAS) whenever R 0 ≤ 1 and that a unique interior endemic equilibrium exists and is GAS if R 0 > 1. We proceed to extend this model with m host species, capturing a class of zoonoses where the cross-species bridge is an arthropod vector. The basic reproduction number of the multi-host-vector, R 0 (m), is derived and shown to be the sum of basic reproduction numbers of the model when each host is isolated with an arthropod vector. It is shown that the disease will persist in all hosts as long as it persists in one host. Moreover, the overall basic reproduction number increases with respect to the host and that bringing the basic reproduction number of each isolated host below unity in each host is not sufficient to eradicate the disease in all hosts. This is a type of "amplification effect," that is, for the considered vector-borne zoonoses, the increase in host diversity increases the basic reproduction number and therefore the disease burden.

Introduction

Zoonoses, infectious diseases caused by pathogens transmissible under natural conditions from vertebrate animals to humans, account for 75% of emerging infectious diseases [START_REF] Taylor | Risk factors for human disease emergence[END_REF], with wildlife being an important source. Zoonotic pathogens represent 61% of all know human pathogens [START_REF] Taylor | Risk factors for human disease emergence[END_REF]. These zoonoses are responsible for over one billion cases of human morbidity, millions of human mortality each year, and over $80 billion USD of global economic burden, including public and animal health and livestock, from 1997 to 2009 [START_REF]World Bank, People, pathogens and our planet, volume 1: Towards a one health approach for controlling zoonotic diseases[END_REF] across the world. Therefore, understanding the dynamics of the underlying mechanisms that drive the inter-host, between hosts, and/or vectors interactions to formulate better control strategies should be a worldwide priority.

Zoonoses are transmitted through four majors transmission routes [START_REF]Zoonotic diseases[END_REF]: i) vector-borne: class that includes diseases such as Lyme disease, tick-borne relapsing fever (TBRF), West Nile virus (WNV), Chagas disease, etc. They are transmitted by the bite of arthropod vectors;

ii) direct transmission: including Brucellosis, hantavirus, influenza, rabies, etc., these zoonoses are due to a close contact with bodily fluids of an infected host, including contacts with fomites;

iii) indirect transmission (e.g., Anthrax, Echinococcosis, Leptospirosis, etc.) that are transmitted through the air by droplet transfer from an infected host to susceptible host; and iv) food-borne or oral transmission: class that includes Toxoplasmosis, Trichinellosis, Salmonellosis, etc. that are caused by ingesting food or water contaminated with a pathogen.

However, according to [START_REF] Johnson | Spillover and pandemic properties of zoonotic viruses with high host plasticity[END_REF], arthropod vectors transmit 40% of zoonoses involving wild animals and 20% of zoonoses involving domestic animals, making vector-borne zoonoses (VBZ) arguably the most important class of zoonoses (see [START_REF] Lloyd-Smith | Epidemic dynamics at the human-animal interface[END_REF] for an excellent review on VBZ). The understanding and control of zoonoses have been hampered because of the complexities of interactions of zoonoses at the interface of humans, animals, and the ecosystem. Mathematical models have long been used to gain insights of key components of the disease in consideration, although most have dealt with only one host and a pathogen [START_REF] Anderson | Infectious Diseases of Humans. Dynamics and Control[END_REF][START_REF] Kermack | Contributions to the mathematical theory of epidemics-i[END_REF] or along with a vector [START_REF] Ross | The prevention of malaria[END_REF][START_REF] Ross | An application of the theory of probabilities to the study of a priori pathometry. part i[END_REF][START_REF] Ross | An application of the theory of probabilities to the study of a priori pathometry. part ii[END_REF][START_REF] Ross | An application of the theory of probabilities to the study of a priori pathometry. part iii[END_REF].

While dynamical models capturing directly transmitted zoonoses, including food-borne and freeliving pathogens, have received a substantial, yet insufficient, attention [START_REF] Begon | Host-host-pathogen models and microbial pest control: the effect of host self regulation[END_REF][START_REF] Begon | Disease and community structure: the importance of host self-regulation in a host-host-pathogen model[END_REF][START_REF] Bowers | A host-host-pathogen model with free-living infective stages, applicable to microbial pest control[END_REF][START_REF] Greenman | Parasite-mediated and direct competition in a two-host shared macroparasite system[END_REF][START_REF] Holt | Infectious disease and species coexistence: a model of lotka-volterra form[END_REF][START_REF] Lloyd-Smith | Epidemic dynamics at the human-animal interface[END_REF] (and the references therein), the literature on the dynamics of multi-host zoonoses is sparse when transmission involves an arthropod vector (mosquitoes, flies, fleas, ticks, etc.) [START_REF] Bowman | A mathematical model for assessing control strategies against West Nile virus[END_REF][START_REF] Cruz-Pacheco | Control measures for chagas disease[END_REF][START_REF]Multi-species interactions in west nile virus infection[END_REF][START_REF] Johnson | Modeling relapsing disease dynamics in a host-vector community[END_REF][START_REF] Simpson | Vector host-feeding preferences drive transmission of multi-host pathogens: West nile virus as a model system[END_REF], most of which are focused on WNV and with at most two hosts. Indeed, even among papers that claim to model a zoonosis, authors usually disregard the multi-host component.

The main goal of this paper is to formulate and investigate the global asymptotic behavior of systems describing the interactions between multiple hosts and vectors that capture the evolution of a class of vector-borne zoonoses. The outline of this endeavor is based on the following steps:

• We begin by the formulating a single host-vector model with an arbitrary number of stages (or relapses) during the host's infectiosity. The proposed model revisits the models in [START_REF] Johnson | Modeling relapsing disease dynamics in a host-vector community[END_REF][START_REF] Palmer | The dynamics of vector-borne relapsing diseases[END_REF] and generalizes them in order to incorporate the heterogeneity of the vector's infection with respect to the host's stage of infection (Section 1).

• We then derive the basic reproduction number and provide a complete investigation of the global stability of equilibria, which surprisingly has not been done (subsections 1.2 and 1.3).

• We formulate a multi-host-vector model where each host's infection dynamics have an arbitrary number of stages (or relapses). The proposed model allows also that some of the hosts are dead-end hosts, that is, they carry the pathogen and may be deadly but do not spread the infection. This feature of the model allows for inspecting whether host diversity increases or mitigates the overall infection.

• We complete the global stability of equilibria of the model with m host species on which a common arthropod vector is feeding. These results are new (section 2).

The basic reproduction number of the multi-host-vector system is R 2 0 (m) = m j=1 R 2 0,j , where R 0,j , for j = 1, 2, . . . , m, is the basic reproduction number of a single host-vector system, when taken in isolation. We prove that the disease-free equilibrium (DFE) is globally asymptotically stable if R 0 (m) ≤ 1. That is, in this case, the disease dies out in all hosts and in the vector populations. If R 0 (m) > 1, a unique interior endemic equilibrium exists and is shown to be globally asymptotically stable. Given the expression of R 2 0 (m) in terms of R 2 0,j , for j = 1, 2, . . . , m, this implies that controlling vector-borne zoonoses requires a coordinated effort in controlling the disease in all hosts. This result echoes the One Health concept [START_REF] Health | [END_REF], which advocates for optimal health for all species and the environment [START_REF]One health, One Health Commission[END_REF]. Moreover, the basic reproduction number R 0 (m) increases with respect to the number of hosts, although the increase of dead-hosts does not affect R 0 (m). This is a type of "amplification effect."

Single host-vector with stage progression

In this section, we derive an SEIR-SI vector-borne disease for a single host-vector and study its asymptotic properties. We later extend to the multi-host case in section 2.

Formulation of the Model

We consider a host population, denoted by N h , that interacts with a population of vectors, denoted by N v . This host-vector interaction leads to a transmission of pathogens from host to vector and conversely. We suppose that the host population is subdivided into susceptible host, denoted by S h , latent host (infected but not infectious) E h , hosts in a successive chain of infected classes I i , for i = 1, 2, . . . , n that represents the different stages of the disease progression, and a class of recovered hosts, denoted by R h . The total host population is therefore,

N h = S h + E h + n i=1 I i + R h . A
host's infection thus follows an SEIR type of framework. Vectors' population is composed of susceptible and infected vectors, denoted respectively by S v and I v , and their infection follows an SI model type of structure. Host's susceptible population is generated via a constant recruitment Λ h and reduced by a per-capita mortality rate of µ h and through an infectious vector bite. The dynamics of the susceptible host is therefore given by:

Ṡh = Λ h -b(N h , N v )β vh S h I v N v -µ h S h , (1.1)
where b(N h , N v ) is the number of mosquito bites per human per unit of time. The new infected hosts become latent and leave this stage either by natural death or entering the infectious stage after an incubation period of 1/ν h . Hence, the dynamics of latent host is given by:

Ėh = b(N h , N v )β vh S h I v N v -(µ h + ν h )E h . (1.2)
The host's incubation period is important for vector-borne diseases in general and vector-borne zoonoses in particular. For instance, the latency period for the Japanese Encephalitis Virus is 10 days in pigs [START_REF] Khan | Dynamics of japanese encephalitis virus transmission among pigs in northwest bangladesh and the potential impact of pig vaccination[END_REF] and 5-15 days [START_REF]Transmission of japanese encephalitis virus, Centers for Disease Control and Prevention[END_REF] for humans. Humans are dead-end host, that is, although deadly, the parasitemia is insufficient to infect biting mosquitoes. Similarly, the incubation period of TBRF is seven days [START_REF] Dworkin | Tick-borne relapsing fever[END_REF][START_REF] Southern | Relapsing fever: a clinical and microbiological review[END_REF]. For other tick-borne diseases, the incubation period varies from more than two weeks and up to nine weeks (see [START_REF]Tickborne diseases of the united states, Centers for Disease Control and Prevention[END_REF] for an overview of the incubation period of tick-borne diseases). The long incubation period of some zoonoses (for instance, in Chagas's disease case, human host may stay asymptomatic his/her whole life), compared to their infectiosity period, make it indispensable to incorporate a latent class in the dynamics of zoonoses.

After the incubation period, the latent host enters the first stage of infectiousness, namely I 1 .

The evolution of the latter follows,

İ1 = ν h E h -(µ h + η 1 + γ 1 )I 1 .
The infected host I 1 leave this stage either by natural death at a rate µ h , through recovery at a per capita rate η 1 , or by progressing to the second stage of infectiousness at a per-capita rate of

γ 1 . Hence, İ2 = γ 1 I 1 -(µ h + η 2 + γ 2 )I 2 .
Similarly, for i = 3, 4, . . . , n -1, the dynamics of the infected at stage i is captured by:

İi = γ i-1 I i-1 -(µ h + η i + γ i )I i , and İn = γ n-1 I n-1 -(µ h + η n )I n .
Hosts' stage-structure infectiousness in vector-borne zoonoses have long been recognized in literature [START_REF] Cruz-Pacheco | Control measures for chagas disease[END_REF]24,[START_REF] Velascohernandez | A model for chagas disease involving transmission by vectors and blood transfusion[END_REF]. Although a review of the epidemiology of VBZ is beyond the scope of this paper, we present in the following a list of some known VB zoonoses to justify the rationale behind incorporating an age-structure host infectiosity in the host's dynamics.

• Tick-borne relapsing fever (TBRF): Caused by Borrelia spirochetes and transmitted by ticks of genus Ornithodoros, the infection has two main stages after the latency period: the febrile and afebrile stages that last on average 3.1 days and 9.25 days, respectively [START_REF] Southern | Relapsing fever: a clinical and microbiological review[END_REF]. Moreover, according to [START_REF] Southern | Relapsing fever: a clinical and microbiological review[END_REF], although there is an average of 3 stages (or relapses), up to 13 relapses could be observed in TBRF. It is worthwhile to note that different nomenclatures have been used to describe different infectious stages during an infection. The terminology of stage progression is used in mathematical modeling [START_REF] Guo | Global dynamics of a staged progression model for infectious diseases[END_REF][START_REF] Guo | Global dynamics of a general class of multistage models for infectious diseases[END_REF][START_REF] Iggidr | multi-compartment models[END_REF], whereas the relapses is used more in ecology [START_REF] Johnson | Modeling relapsing disease dynamics in a host-vector community[END_REF][START_REF] Palmer | The dynamics of vector-borne relapsing diseases[END_REF][START_REF] Southern | Relapsing fever: a clinical and microbiological review[END_REF].

• American Trypanosomiasis: After an incubation period of 1-2 weeks, the infection, also known as Chagas disease (caused by the parasite Trypanosoma cruzi and transmitted by the vector kissing bug), has three main infectious stages: acute (lasts 4-8 weeks), indeterminate, and chronic. In the first two stages, symptoms might be absent, mild or unspecific. Similarly, for the African Trypanosomiasis, also known as sleeping sickness, there are two main infectious stages [START_REF] Franco | Epidemiology of human african trypanosomiasis[END_REF] after a seven year incubation period [START_REF] Wengert | Human african trypanosomiasis with 7-year incubation period: clinical, laboratory and neuroimaging findings[END_REF]. In the first stage, the parasitamia is in the peripheral bloodstream, whereas in the later stage, the parasites enter the nervous system. However, during this long latency period, the parasiteamia in the bloodstream is considered to be sufficient to maintain the transmission cycle [START_REF] Wengert | Human african trypanosomiasis with 7-year incubation period: clinical, laboratory and neuroimaging findings[END_REF].

• West Nile Virus: Transmitted to different species mainly by Culex mosquitoes, WNH has different pathogenicity depending on the host species under consideration. For humans, there is an asymptomatic stage of 3-14 days and two main infectious stages, namely West Nile Fever, and West Nile neuroinvasive disease (WNND).

The recovered hosts are replenished by the sum of the infected of all stages who recovered from the disease and are reduced by natural death rate. As the recovery rate of infectious hosts at stage i is η i , the dynamics of the recovered are given by:

Ṙ = n i=1 η i I i -µ h R.
The susceptible vectors are recruited via a constant recruitment Λ v and die at natural death rate of µ v or via control measures, such as Dichlorodiphenyltrichloroethane (DDT) (for mosquitoes, sandflies, etc.) and Acaricide (for ticks and mites) at a rate of δ v . The infection of susceptible vectors occurs after biting an infectious host. The total number of infectious hosts is

n i=1 I i .
However, for many vector-borne zoonoses, it is reasonable to assume that the vector's infectiousness to hosts' early stages of the infection differs from that of later stages as the parasite load in infectious individuals increases with respect to the duration of infection. Therefore, by denoting β i as the infectiousness of the vector to host of stage i per bite and a as the biting rate, we incorporate the differential infectiousness of vectors to the "age"-structured infection of the host. Hence, the dynamics of the susceptible and infected vectors are respectively given by:

Ṡv = Λ v -aS v n i=1 β i I i N h -(µ v + δ v )S v ,
and

İv = aS v n i=1 β i I i N h -(µ v + δ v )I v .
By the conservation law, the total number of bites on the host by mosquitoes (aN v ) should equal the total number of bites received by host (b

(N h , N v )N h ). Therefore, b(N h , N v )N h = aN v , or equivalently, b(N h , N v ) = aN v N h . (1.3)
By plugging (1.3) into (1.1) and (1.2), the overall dynamics of the host-vector interaction are given by the following system:

                                                               Ṡh = Λ h -a β vh S h I v N h -µ h S h Ėh = a β vh S h I v N h -(µ h + ν h )E h İ1 = ν h E h -(µ h + η 1 + γ 1 )I 1 İ2 = γ 1 I 1 -(µ h + η 2 + γ 2 )I 2 . . . İn-1 = γ n-2 I n-2 -(µ h + η n-1 + γ n-1 )I n-1 İn = γ n-1 I n-1 -(µ h + η n )I n Ṙ = n i=1 η i I i -µ h R Ṡv = Λ v -aS v n i=1 β i I i N h -(µ v + δ v )S v İv = aS v n i=1 β i I i N h -(µ v + δ v )I v (1.4)
The flow diagram of System (1.4) is represented in Fig. [START_REF] Anderson | Infectious Diseases of Humans. Dynamics and Control[END_REF], and the parameters are described in Table 1.1. The subsystem describing the dynamics of the host is triangular, and hence we can disregard the dynamics of the recovered host R. Moreover, the dynamics of the total host and vector populations are given by Ṅh = Λ h -µ h N h and Ṅv = Λ v -(µ h + δ v )N h . Hence, we can deduce that:

lim t→∞ N h = Λ h µ h := Nh and lim t→∞ N v = Λ v µ v + δ v := Nv .
By using the theory of asymptotically autonomous systems for triangular systems [START_REF] Castillo-Chavez | Asymptotically autonomous epidemic models[END_REF][START_REF] Vidyasagar | Decomposition techniques for large-scale systems with nonadditive interactions: Stability and stabilizability[END_REF], we can see that System (1.4) is asymptotically equivalent to:

Host S h E h I 1 I 2 . . . I n R h aβ vh S h I v Nh ν h E h γ 1 I 1 η 1 I 1 η 2 I 2 γ 2 I 2 γ n-1 I n-1 η n I n Λ h µ h S h µ h E h µ h I 1 µ h I 2 µ h I n µ h R h Vectors S v I v Λ v aS v n i=1 β i I i Nh (µ v + δ v )S v (µ v + δ v )I v
                                         Ṡh = Λ h -a β vh S h I v Nh -µ h S h Ėh = a β vh S h I v Nh -(µ h + ν h )E h İ1 = ν h E h -(µ h + η 1 + γ 1 )I 1 İ2 = γ 1 I 1 -(µ h + η 2 + γ 2 )I 2 . . . İn-1 = γ n-2 I n-2 -(µ h + η n-1 + γ n-1 )I n-1 İn = γ n-1 I n-1 -(µ h + η n )I n İv = a( Nv -I v ) n i=1 β i I i Nh -(µ v + δ v )I v (1.5)
Variations of this model have been considered in the literature. Indeed in modeling tick-borne relapsing fever (TBRF), Johnson et al. in [START_REF] Johnson | Modeling relapsing disease dynamics in a host-vector community[END_REF] and Palmer et al. [START_REF] Palmer | The dynamics of vector-borne relapsing diseases[END_REF] considered an SIR-SI model where the dynamics of the infected host is structured into n stages. However, these authors assume that the transmission from vector to host is homogeneous with respect to the host's stage of infection. Authors in [START_REF] Guo | Global dynamics of a staged progression model for infectious diseases[END_REF][START_REF] Guo | Global dynamics of a general class of multistage models for infectious diseases[END_REF] have investigated SIR and SEIR multi-stages models, respectively, although their systems describe infections that are directly transmitted, that is, not coupled with an arthropod vector. Model (1.5) is an SEIR-SI type with heterogeneous infectiousness of vectors. When 1 ν → 0, that is, when the mean latency period is 0, Model (1.5)

is the limit of an SIR-SI model proposed in [START_REF] Johnson | Modeling relapsing disease dynamics in a host-vector community[END_REF][START_REF] Palmer | The dynamics of vector-borne relapsing diseases[END_REF]. Also, Park [START_REF] Park | Age-dependence in Epidemic Models of Vector-borne Infections[END_REF] have proposed an SEIR-SEI stage progression model. In many vector-borne diseases, zoonotic or otherwise, the incubation period of the vector is relatively short compared to the length of infection, as most arthropod vectors die while infected. Therefore, it is reasonable to assume that the dynamics of the vector is captured by an SI structure. In the above mentioned papers, the authors derived the explicit expression of the equilibria and studied their local stability. However, the global analysis of the asymptotic behavior of these models is lacking. In [START_REF] Johnson | Modeling relapsing disease dynamics in a host-vector community[END_REF], the authors assume different natural death rates in each of the infected classes but assume that µ i = 0.01 i + 2 for all i. In Model (1.4), we assumed that all death rates are equal.

However, it is not difficult to transform System (1.4) into a system with different death rates in each classes. Indeed, by choosing

Λ h = µ S S h + µ E E h + n i=1 µ h,i I i + µ R R
, the total population will be constant, and by converting the variables into proportions, we obtain exactly System (1.5).

System (1.5) could be written in a more compact form as follows:

                 Ṡh = Λ h e 1 -a β vh S h I v Nh -µ h S h İh = a β vh S h I v Nh e 1 + A I h İv = a Nv -I v N h 0 β | I h -(µ v + δ v )I v (1.6)
where

I h = (E h , I 1 , . . . , I n ) T , e 1 = (1, 0, . . . , 0) T ∈ IR n+1 , β = (β 1 , β 2 , . . . , β n ) T , x | y is the canonical scalar product (here in IR n+1
), and A is given by:

A =           -(ν + µ h ) 0 . . . . . . 0 ν -α 1 . . . . . . 0 0 γ 1 -α 2 . . . 0 . . . . . . . . . . . . . . . 0 . . . . . . γ n-1 -α n          
,

where α i = γ i + η i + µ h , for i = 1, 2, . . . , n -1 and α n = η n + µ h .
In the following lemma, we prove that the solutions of Model (1.6) remain positive and bounded, that is, the solutions are biologically substantiated.

Lemma 1.1. The region defined by

Ω = (S h , E h , I i , I v ) ∈ IR n+3 + | S h + E h + n i=1 I i ≤ Nh , I v ≤ Nv
is a compact attracting positively invariant set for system (1.6).

The disease-free equilibrium is E 0 = ( Sh , 0 IR n+2 ), where Sh = Λ µ h . This equilibrium always exists and it belongs to Ω.

The infected compartments of System (1.6) are the sum of F and V where

F(I h , I v ) =      a β vh S h I v Nh e 1 a Nv -I v N h 0 β | I h      and V(I h , I v ) = A I h -(µ v + δ v )I v .
Let F = DF(I h , I v ) and V = DV(I h , I v ) be function-valued matrices evaluated at the DFE. We obtain:

V = A 0 0 1,n+1 -(µ v + δ v ) and F =   0 n+1,n+1 aβ vh e 1 aβ hv Nv Nh (0, β) T 0   .
Given the particular form of A, it is not difficult to verify that:

-A -1 =               1 ν + µ h 0 . . . . . . 0 ν α 1 (ν + µ h ) 1 α 1 . . . . . . 0 γ 1 ν α 1 α 2 (ν + µ h ) γ 1 α 1 α 2 1 α 2 . . . 0 . . . . . . . . . . . . . . . γ 1 . . . γ n-1 ν α 1 . . . α n (ν + µ h ) • • • • • • γ n-1 α n-1 α n 1 α n               (1.7)
Hence, the next generation matrix is given by:

-F V -1 =    0 n+1,n+1 a β vh µ v + δ e 1 -a Nv Nh (0, β) T A -1 0    .
The basic reproduction number, defined as the average number of secondary cases produced by an infected host during its infectious period while interacting with a completely susceptible population, is computed using the next generation method [START_REF] Diekmann | Mathematical epidemiology of infectious diseases[END_REF][START_REF] Van Den Driessche | reproduction numbers and sub-threshold endemic equilibria for compartmental models of disease transmission[END_REF]. Denoted by R 0 , it is the spectral radius of the next generation matrix -F V -1 . Therefore, for Model (1.6), the basic reproduction number is,

R 2 0 = a 2 β vh µ v + δ Nv Nh (0, β) T (-A -1 )e 1 := a 2 β vh Nv ν Nh (µ v + δ)(ν + µ h ) n i=1 γ 1 γ 2 . . . γ i-1 α 1 α 2 . . . α i β i with γ 0 = 1.
This quantity is intrinsically tied to the dynamics of the disease as we will discuss in the next subsection.

Global stability of the DFE

By the derivation of the basic reproduction number, the DFE is locally asymptotically stable if R 0 < 1 and unstable if R 0 > 1 [START_REF] Diekmann | Mathematical epidemiology of infectious diseases[END_REF][START_REF] Van Den Driessche | reproduction numbers and sub-threshold endemic equilibria for compartmental models of disease transmission[END_REF]. The following theorem provides a global result.

Theorem 1.1.

The DFE is globally asymptotically stable if R 2 0 ≤ 1.

Proof.

The instability of the DFE when R 0 > 1 follows from [START_REF] Van Den Driessche | reproduction numbers and sub-threshold endemic equilibria for compartmental models of disease transmission[END_REF]. As for the global stability of the DFE when R 0 ≤ 1, let us consider the candidate Lyapunov function:

V(I h , I v ) = L(I h ) + c v I v where L(I h ) = c 0 E h + n i=1 c i I i , with c v = Nh a Nv
and the coefficients c i for i = 1, 2, . . . , n -1 are given by the following induction relationship:

                   c 0 = ν ν + µ h c 1 c 1 = β 1 α 1 + γ 1 α 1 α 2 β 2 + • • • + γ 1 . . . γ n-1 α 1 α 2 . . . α n β n , c i+1 = 1 γ i (α i c i -β i ) , for i = 1, 2, . . . , n -2, c n = β n α n .
(

The derivative along the trajectories of the system is:

V(I h , I v ) = L(I h ) + c v İv = c 0 Ėh + n i=1 c i İi + c v İv = c 0 a β vh S h I v Nh -(µ h + ν h )E h + c 1 (ν h E h -α 1 I 1 ) + n i=2 c i İi + c v İv = c 0 a β vh S h I v Nh -c 1 α 1 I 1 + c 2 (γ 1 I 1 -α 2 I 2 ) + n i=3 c i İi + c v İv = c 0 a β vh S h I v Nh -c 1 α 1 I 1 + (c 1 α 1 -β 1 )I 1 -c 2 α 2 I 2 + n i=3 c i İi + c v İv , since c 2 γ 1 = α 1 c 1 -1. Therefore, V(I h , I v ) = c 0 a β vh S h I v Nh -β 1 I 1 -c 2 α 2 I 2 + n i=3 c i İi + c v İv .
Similarly, by using successively the induction relationship (1.8), we obtain,

V(I h , I v ) = c 0 a β vh S h I v Nh - n i=1 β i I i + c v a(N v -I v ) n i=1 β i I i Nh -(µ v + δ v )I v . (1.9)
Considering the expression of c v = Nh a Nv and the inequality S h ≤ Nh , Equality (1.9) leads to:

V(I h , I v ) = c 0 aβ vh S h Nh I v -I v n i=1 β i I i - Nh (µ v + δ v ) a Nv I v = c 0 aβ vh S h Nh - Nh (µ v + δ v ) a Nv I v -I v n i=1 β i I i = Nh (µ v + δ v ) a Nv a 2 β vh Nv c 0 (µ v + δ v ) Nh S h Nh -1 I v -I v n i=1 β i I i = Nh (µ v + δ v ) a Nv R 2 0 S h Nh -1 I v -I v n i=1 β i I i , since c 0 = ν ν + µ h n i=1 γ 1 γ 2 . . . γ i-1 α 1 α 2 . . . α i β i . Thus, we conclude that V(I h , I v ) ≤ 0, whenever R 2 0 ≤ 1 as S h ≤ Nh . Moreover, V(I h , I v ) = 0 if I v = 0 or if R 2 0 = 1, S h = Nh and n i=1
β i I i = 0. Hence, we deduce that the largest set contained in

{(I h , I v ) ∈ Ω | V(I h , I v ) = 0} is
reduced to the DFE. Since the set Ω is compact and positively invariant, by LaSalle's invariance principle [START_REF] Bhatia | Stability Theory of Dynamical Systems[END_REF][START_REF] Lasalle | Stability by Liapunov's direct method[END_REF], the DFE is globally asymptotically stable on Ω. The attractiveness of Ω makes the DFE globally asymptotically stable on IR n+3 .

Endemic Equilibria

In this section, we discuss the existence conditions for the endemic equilibrium and investigate its global stability whenever it exists. Let (S * h , I * h , I * v ) T 0 be an endemic equilibrium for Model (1.6). Therefore (S * h , I * h , I * v ) T satisfies the endemic relations:

                                                     0 = Λ h -a β vh S h * I * v Nh -µ h S * h 0 = a β vh S h I * v Nh -(µ h + ν h )E * h E * h = α 1 ν h I * 1 I * 2 = γ 1 α 2 I * 1 I * 3 = γ 2 α 3 I * 2 = γ 1 γ 2 α 2 α 3 I * 1 . . . I * n-1 = γ n-2 α n-1 I * n-2 = γ 1 γ 2 . . . γ n-2 α 2 α 3 . . . α n-1 I * 1 I * n = γ n-1 α n I * n-1 = γ 1 γ 2 . . . γ n-1 α 2 α 3 . . . α n I * 1 0 = a( Nv -I * v ) n i=1 β i I * i Nh -(µ v + δ v )I * v (1.10)
where α i = µ i + η i + γ i . It follows from (1.10) that there is an interior endemic equilibrium if I * 1 > 0. Moreover, we can deduce that:

n i=1 β i I * i = I * 1 β 1 + γ 1 α 2 β 2 + • • • + γ 1 γ 2 . . . γ n-1 α 2 . . . α n β n = α 1 c 1 I * 1 . (1.11)
The first equation of (1.10) implies that:

S * h = µ h Nh µ h + aβ vh I * v Nh and E * h = aβ vh µ h + ν h µ h I * v µ h + aβ vh I * v Nh .
We have

I * 1 = µ h (µ v + δ v )(R 2 0 -1) aµ h Nh ξ + a 2 β vh Nv ξ , (1.12)
where The endemic equilibrium for the multi-host Model (1.6) is globally asymptotically stable whenever R 2 0 > 1.

ξ = β 1 + γ 1 c 2 = α 1 c 1 ,

Proof.

Let us consider the Lyapunov function

W = b 0 (S h -S * h log S h +E h -E * h log E h )+ n i=1 b i (I i -I * i log I i )+v 1 (S v -S * v log S v +I v -I * v log I v )+K
where

K = v 1 (S * v -S * v log S * v + I * v -I * v log I * v ) + b 0 (S * h -S * h log S * h + E * h -E * h log E * h ) + n i=1 b i (I * i -I * i log I * i )
and b i , for i = 0, 1, . . . , n, and v 1 are positive constants. The function W is therefore positivedefinite with respect to the endemic equilibrium. The challenge is how to choose these coefficients to make the derivative of W along the trajectories of (1.6) negative-definite. To ease the notations and WLOG, we denote βi = aβ i 1 Nh , βh = ( β1 , . . . , βn ) T and βv = aβ vh 1 Nh .

The derivation of W along the trajectories of (1.6) is as follows:

Ẇ = b 0 1 - S * h S h Λ h -βv S h I v -µ h S h + b 0 1 - E * h E h βv S h I v -(µ h + ν h )E h +b 1 1 - I * 1 I 1 (ν h E h -(µ h + η 1 + γ 1 )I 1 ) + n i=2 b i 1 - I * i I i (γ i-1 I i-1 -(µ h + η i + γ i )I i ) +v 1 1 - S * v S v Λ v -S v n i=1 βi I i -(µ v + δ v )S v +v 1 1 - I * v I v S v n i=1 βi I i -(µ v + δ v )I v .
(1.13) By grouping all linear terms of (1.13) in E h , I i , for i = 1, 2, . . . , n and I v , and denoting it by L, we obtain after lengthy computations and arrangements:

L = b T Ã + v 1 S * v βT h I h + (b 1 ν h -b 0 (µ h + ν h ))E h + b 0 βv S * h -v 1 (µ v + δ v ) I v (1.
= -v 1 S * v Ã-T βh and b 0 = ν h µ h + ν h b 1 .
More precisely, from the expression of b, and using endemic relations (1.10), the first component of b could take the following different forms:

b 1 = v 1 S * v n i=1 γ 1 ...γ i-1 α 1 ...α i βi = v 1 S * v 1 α 1 β1 + γ 1 α 1 α 2 β2 + • • • + γ 1 . . . γ n-1 α 1 α 2 . . . α n βn = v 1 S * v 1 α 1 β1 + γ 1 α 2 β2 + • • • + γ 1 . . . γ n-1 α 2 . . . α n βn = v 1 S * v 1 α 1 n i=1 βi I * i I * 1 . (1.15)
Also, b 0 could clearly be expressed as

b 0 = ν h µ h + ν h b 1 = ν h µ h + ν h v 1 S * v 1 α 1 n i=1 βi I * i I * 1 .
(1.16)

The use of (1.16) cancels the coefficient of

I v in (1.14). Indeed, b 0 βv S * h -v 1 (µ v + δ v ) = ν h µ h + ν h v 1 S * v 1 α 1 n i=1 βi I * i I * 1 βv S * h -v 1 (µ v + δ v ) = v 1 ν h (µ h + ν h )α 1 I * 1 S * v n i=1 βi I * i βv S * h -v 1 (µ v + δ v ) = v 1 ν h (µ h + ν h )α 1 I * 1 (µ v + δ v ) I * v βv S * h -v 1 (µ v + δ v ) = v 1 ν h (µ h + ν h )α 1 I * 1 (µ v + δ v ) (µ h + ν h )E * h -v 1 (µ v + δ v ) = v 1 ν h (µ h + ν h )α 1 I * 1 (µ v + δ v ) (µ h + ν h ) α 1 I * 1 ν h -v 1 (µ v + δ v ) = 0.
Therefore, with this choice for b and b 0 , all linear terms in (1.13) cancel. The remaining terms in Ẇ are:

Ẇ = v 1 1 - S * v S v (Λ v -(µ v + δ v )S v ) + v 1 (µ v + δ v )I * v -v 1 I * v I v S v n i=1 βi I i +b 0 1 - S * h S h (Λ h -µ h S h ) + b 0 (µ h + ν h )E * h -βv S h I v E * h E h + b 1 (µ h + η 1 + γ 1 )I * 1 -b 1 ν h I * 1 E h I 1 + n-1 i=2 b i (µ h + η i + γ i )I * i -b i γ i-1 I * i I i-1 I i + b n (µ h + η n )I * n -b n γ n-1 I * n I n-1 I n = v 1 (µ v + δ v )S * v 2 - S * v S v - S v S * v Av +2v 1 (µ v + δ v )I * v -v 1 (µ v + δ v ) I * v S * v S v -v 1 I * v I v S v n i=1 βi I i +b 0 µ h S * h 2 - S * h S h - S h S * h A h +2b 0 (µ h + ν h )E * h -b 0 (µ h + ν h ) E * h S * h S h + b 0 -βv S h I v E * h E h +b 1 (µ h + η 1 + γ 1 )I * 1 -b 1 ν h I * 1 E h I 1 + n-1 i=2 b i (µ h + η i + γ i )I * i -b i γ i-1 I * i I i-1 I i +b n (µ h + η n )I * n -b n γ n-1 I * n I n-1 I n .
By using the endemic relations and (1.15-1.16), b 0 and b 1 could be written as:

b 0 = ν h µ h + ν h b 1 = ν h µ h + ν h v 1 S * v 1 α 1 n i=1 βi I * i I * 1 = v 1 µ h + ν h S * v n i=1 βi I * i E * h = v 1 βv S * v n i=1 βi I * i S * h I * v .
Hence,

Ẇ = v 1 A v + 2v 1 (µ v + δ v )I * v -v 1 (µ v + δ v ) I * v S * v S v -v 1 I * v I v S v n i=1 βi I i + b 0 A h + 2b 0 (µ h + ν h )E * h -b 0 (µ h + ν h ) E * h S * h S h -b 0 βv S h E * h I v E h +b 1 (µ h + η 1 + γ 1 )I * 1 -b 1 ν h I * 1 E h I 1 + n-1 i=2 b i (µ h + η i + γ i )I * i -b i γ i-1 I * i I i-1 I i +b n (µ h + η n )I * n -b n γ n-1 I * n I n-1 I n .
Taking into account the construction of b i , we have

2b 0 (µ h + ν h )E * h = 2v 1 (µ v + δ v )I * v = 2v 1 S * v n i=1 βi I * i , b i = v 1 S * v n j=i βj I * j (µ h + η i + γ i )I * i = v 1 S * v n j=i βj I * j α i I * i , ∀i = 1, 2, . . . , n -1,
and,

b n = v 1 βn S * v . It follows that: Ẇ/v 1 = A v + (b 0 A h /v 1 ) + 2S * v n i=1 βi I * i -S * v n i=1 βi I * i S * v S v - I * v I v S v n i=1 βi I i + 2S * v n i=1 βi I * i -S * v n i=1 βi I * i S * h S h -S * v n i=1 βi I * i S * h I * v S h E * h I v E h + S * v n i=1 βi I * i -S * v n i=1 βi I * i E * h I * 1 E h I 1 + n-1 i=2   S * v n j=i βj I * j -S * v n j=i βj I * j I * i-1 I * i I i-1 I i   + βn S * v I * n -βn S * v I * n I * n-1 I * n I n-1 I n .
By separating the first index from the sum, we obtain:

Ẇ/v 1 = A v + (b 0 A h /v 1 ) + 5S * v β1 I * 1 -S * v β1 I * 1 S * v S v - I * v I v S v β1 I 1 -S * v β1 I * 1 S * h S h -S * v β1 I * 1 S h E * h I v S * h I * v E h -S * v β1 I * 1 I * 1 E h E * h I 1 +5S * v n i=2 βi I * i -S * v n i=2 βi I * i S * v S v - I * v I v S v n i=2 βi I i -S * v n i=2 βi I * i S * h S h -S * v n i=2 βi I * i S h E * h I v S * h I * v E h -S * v n i=2 βi I * i I * 1 E h E * h I 1 + n i=2   S * v n j=i βj I * j -S * v n j=i βj I * j I * i I i-1 I i I * i-1   .
Therefore,

Ẇ/v 1 = A v + (b 0 A h /v 1 ) + S * v β1 I * 1 5 - S * v S v - S v S * v I 1 I * 1 I * v I v - S * h S h - S h I v E * h S * h I * v E h - E h I * 1 E * h I 1 +S * v n i=2 βi I * i 5 - S * v S v - S v I i I * v S * v I * i I v - S * h S h - S h I v E * h S * h I * v E h - E h I * 1 E * h I 1 + n i=2   S * v n j=i βj I * j -S * v n j=i βj I * j I * i I i-1 I i , I * i-1   .

By using the properties of nested sums,

I * i I i-1 I i I * i-1   = n i=2 βi I * i i i=2 I * i I i-1 I i I * i-1 . Finally, Ẇ/v 1 = A v + (b 0 A h /v 1 ) + S * v β1 I * 1 5 - S * v S v - S v S * v I 1 I * 1 I * v I v - S * h S h - S h I v E * h S * h I * v E h - E h I * 1 E * h I 1 +S * v n i=2 βi I * i 4 + i - S * v S v - S v I i I * v S * v I * i I v - S * h S h - S h I v E * h S * h I * v E h - E h I * 1 E * h I 1 - i i=2 I * i I i-1 I i I * i-1
.

By setting v 1 = 1 and replacing A v and A h by their respective values, the final expression of Ẇ is:

Ẇ = (µ v + δ v )S * v 2 - S * v S v - S v S * v + b 0 µ h S * h 2 - S * h S h - S h S * h +S * v β1 I * 1 5 - S * v S v - S v S * v I 1 I * 1 I * v I v - S * h S h - S h I v E * h S * h I * v E h - E h I * 1 E * h I 1 +S * v n i=2 βi I * i 4 + i - S * v S v - S v I i I * v S * v I * i I v - S * h S h - S h I v E * h S * h I * v E h - E h I * 1 E * h I 1 - i i=2 I * i I i-1 I i I * i-1
(1.17)

with b 0 = S * v n i=1 βi I * i (µ h + ν h )E * h
. By using the geometric-arithmetic mean inequality, we conclude that Ẇ is negative definite, which proves the global asymptotic stability of the endemic equilibrium.

This result of global stability of the endemic equilibrium is new and extend the results of [START_REF] Palmer | The dynamics of vector-borne relapsing diseases[END_REF] in which the endemic equilibrium is proven to be locally asymptotically stable if R 2 0 > 1 and R 2 0 is close to one. Theorems 1.1 and 1.2 solve completely the global dynamics of the vector-borne staged-progression model (1.6). The basic reproduction number R 2 0 acts as a sharp threshold for the disease in the sense that the latter dies out if R 2 0 ≤ 1 and persists whenever R 2 0 > 1.

Now that the dynamics of the single host-vector staged progression is solved, we tackle the dynamics of multi-host vector-borne zoonoses in the next section.

Multi-host vector-borne models

Vector-borne zoonoses are central in understanding the dynamics of zoonoses in general as vectors play the bridging role in transporting the infection from one species to another. Many zoonotic pathogens are shared in multiple hosts. For instance, for West Nile virus, the pathogen has been found in more 300 bird species as well as in other mammals such as horses, bats, and squirrels, among others [START_REF] Marfin | Widespread west nile virus activity, eastern united states, 2000[END_REF]. Therefore, it is important to investigate a general model that captures the dynamics of an arbitrary number of host species interacting with an arthropod vector that bridges the infection among hosts. In [START_REF] Cruz-Pacheco | Control measures for chagas disease[END_REF][START_REF]Multi-species interactions in west nile virus infection[END_REF], the authors investigated the dynamics of Chagas' disease and WNV, respectively, with two host species. In modeling TBRF, Johnson et al. [START_REF] Johnson | Modeling relapsing disease dynamics in a host-vector community[END_REF] proposed a two hosts and one vector model, namely pine squirrels and deer mice. These two hosts are bitten by the same vector (ticks).

Building upon the same scheme as in Model (1.4), we derive a multi-host vector borne zoonoses model with m hosts and n stages of infection. However, the infection terms in the multi-host necessitates a detailed formulation as it differs from the one host scenario. For instance, the equation of susceptible of Host j is given by

Ṡj = Λ j -b j (N j , N v )β vh S j I v N v -µ h S h , (2.1)
where b j (N j , N v ) is the number of bites per Host j per unit of time. Hence, the total number of bites on Host j is b j (N j , N v )N j . The equation of susceptible vectors is given by:

Ṡv = Λ v - m j=1 a j β hv S v n i=1 β j,i I j,i N j -(µ v + δ v )S v
where a j is the number of bites per mosquito per unit of time on Host j.

By the conservation law, the total number of bites on Host j by mosquitoes (a j N v ) should equal the total number of bites received by Host j (b j (N j , N v )N j ). Therefore,

b j (N j , N v )N j = a j N v , or equivalently, b j (N j , N v ) = a j N v N j .
Hence, for j = 1, 2, . . . , m, Equation (2.1) leads to,

Ṡj = Λ j -a j β vh S j I v N j -µ j S j .
Therefore, the dynamics of a model that captures the interactions between m host species and an arthropod vector is given by:

                                                             For j = 1, 2 . . . , m : Ṡj = Λ j -a j β v,j S j I v Nj -µ j S j Ėj = a j β v,j S j I v Nj -(µ j + ν j )E h,j İj,1 = ν j E h,j -(µ j + η j,1 + γ j,1 )I j,1 İj,2 = γ j,1 I j,1 -(µ j + η j,2 + γ j,2 )I j,2 . . . İj,n-1 = γ j,n-2 I j,n-2 -(µ j + η j,n-1 + γ j,n-1 )I j,n-1 İj,n = γ j,n-1 I j,n-1 -(µ j + η j,n + γ j,n )I j,n Ṡv = Λ v - m j=1 a j S v n i=1 β j,i I j,i Nj -(µ v + δ v )S v İv = m j=1 a j S v n i=1 β j,i I j,i Nj -(µ v + δ v )I v . (2.2)
Parameters in System (2.2) are described in Table 2. Remark 2.1. Some authors [START_REF] Bowman | A mathematical model for assessing control strategies against West Nile virus[END_REF][START_REF] Dobson | Population dynamics of pathogens with multiple host species[END_REF][START_REF] Simpson | Vector host-feeding preferences drive transmission of multi-host pathogens: West nile virus as a model system[END_REF] assumed that the vectors' bitting rate is constant on all hosts (that is, vectors bite a certain fixed number of bites a day regardless of hosts) and denoted by a. By the conservation law,

b j (N j , N v ) = aN v m j=1 N j .
Our approach assumes that a j represents the bitting rate of vector on Host j. This differential vector biting rates on hosts embodies also the well-documented bitting/feeding preference of vectors with respect to the hosts ( [START_REF] Johnson | Modeling relapsing disease dynamics in a host-vector community[END_REF][START_REF] Simpson | Vector host-feeding preferences drive transmission of multi-host pathogens: West nile virus as a model system[END_REF]). However, the two approaches are mathematically equivalent if the host populations are asymptotically constant.

System (2.2) models a range of modeling scenarios in vector-borne zoonoses. For instance, if a subset K of the considered hosts carries the pathogen but are dead-end hosts, that is, they do not spread the infection, it is sufficient to let β k,i = 0 for all k ∈ K and i = 1, 2, . . . , n. It also captures the case where different hosts have different epidemiological responses to the infection. That is, by appropriately choosing the parameters in System (2.2), the model could describe the case the infection follows an SIR-type of structure for a collection of hosts and SEIR-type or SI-type with multiple stage of infections for other hosts.

With β j = (β j,1 , . . . , β j,m ) T , System (2.2) could be written in a more compact vectorial form as:

                                                           Ṡ1 = Λ 1 e 1 -a 1 β v,1 S 1 I v N1 -µ 1 S 1 Ṡ2 = Λ 2 e 1 -a 2 β v,2 S 2 I v N2 -µ 2 S 2 . . . Ṡm = Λ m e 1 -a m β v,m S m I v Nm -µ m S m İ1 = a 1 β v,1 S 1 I v N1 e 1 + A 1 I 1 İ2 = a 2 β v,2 S 2 I v N2 e 1 + A 2 I 2 . . . İm = a m β v,m S m I v Nm e 1 + A m I m İv = m j=1 a j Nv -I v Nh,j 0
β j | I j -(µ v + δ v )I v , (2.3)
where I j = (E j , I j,1 , I j,2 . . . , I j,n ) T is the vector of infected of Host j (j = 1, 2, . . . , m) and

A j =           -(ν j + µ j ) 0 . . . . . . 0 ν j -α j,1 . . . . . . 0 0 γ j,1 -α j,2 . . . 0 . . . . . . . . . . . . . . . 0 . . . . . . γ j,n-1 -α j,n          
where α j,i = γ j,i + η j,i + µ j , for j = 1, . . . , m and i = 1, 2, . . . , n. The term α j,i represents the mean period for which infected population of Host j in stage i leave this stage.

Basic Reproduction Number and Basic Properties

In this subsection, we derive the basic reproduction number of the multi-host system. First of all, the solutions of Model (2.3) are bounded as a result of the following Lemma, whose proof is straightforward.

Lemma 2.1. The region defined by

Ω = (S j , E j , I j,i , I v ) ∈ IR m(n+2)+1 + | S j + E j + n i=1 I j,i ≤ Λ j µ j , for j = 1, . . . m, I v ≤ Λ v µ v
is a compact attracting positively invariant set for system (2.3).

The disease-free equilibrium of Model (2.3) is (S * , 0

IR (n+1)m+1 ) with S * = [S * 1 , S * 2 , . . . , S * m ] = Λ 1 µ 1 , Λ 2 µ 2 , . . . , Λ m µ m
. The vector field of Model (2.3) could be subdivided into new infections vector F and transition vector V where

F(I 1 , . . . , I m , I v ) =                a 1 β v,1 S 1 I v N1 e 1 a 2 β v,2 S 2 I v N2 e 1 . . . a m β v,m S m I v Nm e 1 m j=1 a j Nv -I v Nh,j 0
β j | I j                , V(I 1 , . . . , I m , I v ) =           A 1 I 1 A 2 I 2 . . . A m I m -(µ v + δ v )I v          
.

The Jacobian matrices at the DFE of F and

V are F = DF(I 1 , I 2 . . . . , I m , I v ) DFE and V = DV(I 1 , I 2 . . . . , I m , I v ) DFE
and are given by:

F =            0 n+1,n+1 0 n+1,n+1 • • • • • • a 1 β v,1 e 1 0 n+1,n+1 0 n+1,n+1 • • • • • • a 2 β v,2 e 1 . . . . . . . . . . . . . . . 0 n+1,n+1 0 n+1,n+1 • • • • • • a m β v,m e 1 a 1 Nv N1 (0, β 1 ) T a 2 Nv N2 (0, β 2 ) T • • • a m Nv Nm (0, β m ) T 0            , and V = diag(A 1 , A 2 , . . . , A m , -µ v -δ v ).
It follows that

-V -1 = diag -A -1 1 , -A -1 2 , . . . , -A -1 m , 1 µ v + δ v .
Thus, the next generation matrix is:

-F V -1 =               0 n+1,n+1 • • • • • • a 1 β v,1 µ v + δ v e 1 0 n+1,n+1 • • • • • • a 2 β v,2 µ v + δ v e 1 . . . . . . . . . . . . 0 n+1,n+1 • • • • • • a m β v,m µ v + δ v e 1 a 1 Nv N1 (0, β 1 ) T (-A -1 1 ) • • • a m Nv Nm (0, β m ) T (-A -1 m ) 0              
.

The matrix (-F V -1 ) 2 is of rank one, and its largest eigenvalue is:

R 2 0 (m) = m j=1 a 2 j β v,j µ v + δ v Nv Nh,j (0, β j ) T (-A -1 j )e 1 := m j=1 R 2 0,j ,
where R 2 0,j is the basic reproduction number when Host j is the only coupled host to the vector. A more explicit expression of R 2 0 could be obtained by computing A -1 j . Indeed,

-A -1 j =               1 ν j + µ j 0 . . . . . . 0 ν j α j,1 (ν j + µ j ) 1 α j,1 . . . . . . 0 γ j,1 ν α j,1 α j,2 (ν j + µ j ) γ j,1 α j,1 α j,2 1 α j,2 . . . 0 . . . . . . . . . . . . . . . γ j,1 . . . γ j,n-1 ν j α j,1 α j,2 . . . α j,n (ν j + µ j ) • • • • • • γ j,n-1 α j,n-1 α j,n 1 α j,n               . ( 2.4)
Hence, we deduce that, (0,

β j ) T (-A -1 j )e 1 = β j,1 ν j α j,1 (ν j + µ j ) + β j,2 γ j,1 ν j α j,1 α j,2 (ν j + µ j ) + . . . +β j,n γ j,1 . . . γ j,n-1 ν j α j,1 α j,2 . . . α j,n (ν j + µ j ) = ν j ν j + µ j β j,1 α j,1 + β j,2 γ j,1 α j,1 α j,2 + • • • + β j,n γ j,1 . . . γ j,n-1 α j,1 α j,2 . . . α j,n = ν j ν j + µ j n i=1 γ j,1 γ j,2 . . . γ j,i-1 α j,1 α j,2 . . . α j,i β j,i with γ j,0 = 1.
Finally, the basic reproduction number of a multi-host-vector model with n stages is given by

R 2 0 (m) = m j=1 R 2 0,j ,
where

R 2 0,j = a 2 j β v,j µ v + δ v Nv Nh,j ν j ν j + µ j n i=1 γ j,1 γ j,2 . . . γ j,i-1 α j,1 α j,2 . . . α j,i β j,i .
This result generalizes the results of [START_REF] Cruz-Pacheco | Control measures for chagas disease[END_REF][START_REF]Multi-species interactions in west nile virus infection[END_REF][START_REF] Johnson | Modeling relapsing disease dynamics in a host-vector community[END_REF], for which j = 2, and [START_REF] Palmer | The dynamics of vector-borne relapsing diseases[END_REF] for which j = 1. Similar remarks hold for [START_REF] Bowman | A mathematical model for assessing control strategies against West Nile virus[END_REF].

The overall basic reproduction number increases with respect to the number of hosts. Moreover, the addition of dead-end hosts, for which the basic reproduction number is zero in isolation, keeps the overall basic reproduction number steady. A similar, yet contrasting, result has been obtained in [START_REF] Dobson | Population dynamics of pathogens with multiple host species[END_REF]. Indeed, in [START_REF] Dobson | Population dynamics of pathogens with multiple host species[END_REF] the authors found that the basic reproduction number in multi-host models, when not coupled with an arthropod vector, increases with host diversity for densitydependent transmission and decreases for frequency-dependent transmission. Our model consists of frequency-dependent transmission, but R 2 0 (m) increases with host diversity. This contrast from [START_REF] Dobson | Population dynamics of pathogens with multiple host species[END_REF] stems from our modeling approach, particularly the total number of vector bites on hosts. Indeed, as pointed out in Remark 2.1, we consider a frequency-dependent infection process but in which the arthropod vector has a fixed number of bites on each Host j, that is b j (N j , N v ). In [START_REF] Bowman | A mathematical model for assessing control strategies against West Nile virus[END_REF][START_REF] Dobson | Population dynamics of pathogens with multiple host species[END_REF][START_REF] Simpson | Vector host-feeding preferences drive transmission of multi-host pathogens: West nile virus as a model system[END_REF], the authors define b( m j=1 N j , N v ) as the number of bites on all hosts. The latter will lead to a factor of m j=1 N j in the denominator of the infection term. Thus increasing host diversity decreases the infection force and therefore the basic reproduction number. In spite of this difference on the overall basic reproduction numbers, the systems are mathematically asymptotically equivalent if the host populations are either constant or asymptotically constant. Moreover, since the hosts populations are asymptotically constant, our transmission term could be seen as density-dependent, for which the monotonicity of R 2 0 (m) coincides with that of [START_REF] Dobson | Population dynamics of pathogens with multiple host species[END_REF].

Global Stability of Equilibria

In this subsection, we study the global behavior of the steady states of models describing VBZ, that is Model (2.3), with respect to the overall basic reproduction number R 2 0 (m).

Theorem 2.1.

The DFE is globally asymptotically stable for the multi-host Model (2.3) whenever R 2 0 (m) ≤ 1.

Proof. We consider a barycentric-type Lyapunov function, a weighted sum of Lyapunov functions of the one-host case:

V = m j=1 p j L j + I v
where L j = c j,0 E j + c j,1 I j,1 + c j,2 I j,2 + • • • + c j,n I j,n and p j = a j β j,v

Nv Nh,j . The coefficients c j,i for i = 1, 2, . . . , n are similar to those chosen in the proof of Theorem 1.1. They are defined recursively as follows:

                     c j,0 = ν j ν j + µ j c j,1 c j,1 = β j,1 α j,1 + γ j,1 α j,1 α j,2 β j,2 + • • • + γ j,1 . . . γ j,n-1 α j,1 α j,2 . . . α j,n β j,n c j,i+1 = 1 γ j,i (α j,i c j,i -β j,i ) , ∀i = 1, . . . , n -2 c j,n = β j,n α j,n .
(2.5)

The derivative of V along the solutions of System (2. By using the induction relationship (2.5) and as in the proof of Theorem 1.1, we obtain

Lj = c j,0 a j β v,j S j I v Nj - n i=1 β j,i I j,i . (2.6)
Thus, the derivative of the function V becomes:

V = m j=1 p j c j,0 a j β v,j S j I v Nj - n i=1 β j,i I j,i +   m j=1 a j ( Nv -I v ) n i=1 β j,i I j,i Nj -(µ v + δ v )I v   = m j=1 a j Nv Nj -p j n i=1 β j,i I j,i - m j=1 a j I v n i=1 β j,i I j,i Nj +   m j=1 p j c j,0 a j β v,j S j Nj -(µ v + δ v )   I v = - m j=1 a j I v n i=1 β j,i I j,i Nj +   m j=1 c j,0 a 2 j β v,j
Nv Nj

S j Nj -(µ v + δ v )   I v , since p j = a j Nv Nj .
Hence, we obtain,

V = -m j=1 a j β j,v I v n i=1 I j,i Nj + (µ v + δ v ) m j=1 R 2 0,j S j Nj -1 I v ≤ -I v m j=1 a j β j,v n i=1 I j,i Nj + (µ v + δ v ) R 2 0 (m) -1 I v , since S j ≤ Nj , ∀j = 1, . . . , m, ≤ 0 whenever R 2 0 (m) ≤ 1.
(2.7)

It is not difficult to see that the largest invariant set within { V = 0} is reduced to the DFE, which is therefore globally asymptotically stable in the compact set Ω thanks to ([4], Theorem 3.7.11, page 346, [START_REF] Lasalle | Stability theory for ordinary differential equations[END_REF] Theorem 3). Since Ω is an attractive set, it follows that the DFE is globally asymptotically stable in IR m(n+2)+1 .

Theorem 2.2.

A unique endemic equilibrium exists for the multi-host Model (2.3) whenever R 2 0 (m) > 1.

Proof.

The endemic relations of (2.3) are as follows, for j = 1, 2, . . . , m:

E j = α j,1 ν j I * j,1 , I * j,1 = α j,2 γ j,1 I * j,2 ,
and for i = 2, . . . , n,

I * j,i-1 = α j,i γ j,i-1 I * j,i .
Hence, as in (1.11) for the one host case, we can show that:

n i=1 β j,i I * j,i = α j,1 c j,1 I * j,1 . (2.8)
Indeed, the relationship (2.8) could be shown by expressing all components of the endemic equilibrium in terms of I * v . From the first equation of (2.2), taken at the equilibrium, we obtain:

S * j = µ j N 2 j µ j Nj + a j β v,j I * v E * j = a j β v,j µ j + ν j µ j N 2 j µ j Nj + a j β v,j I * v I * v Nj = a j β v,j µ j + ν j µ j Nj I * v µ j Nj + a j β v,j I * v I * j,1 = ν j α j,1 E * j = ν j α j,1 a j β v,j µ j + ν j µ j Nj I * v µ j Nj + a j β v,j I * v .
We can express, similarly, the expressions of I * j,i in term of I * v and by summing them up, we obtain (2.8).

From the equation of vectors in (2.2) at the equilibrium, we obtain

0 = ( Nv -I * v ) m j=1 a j n i=1 β j,i I * j,i Nj -(µ v + δ v )I * v = ( Nv -I * v ) m j=1 a j α j,1 c j,1 I * j,1 Nj -(µ v + δ v )I * v = ( Nv -I * v ) m j=1 a j α j,1 c j,1 Nj ν j α j,1 a j β v,j µ j + ν j µ j Nj I * v µ j Nj + a j β v,j I * v -(µ v + δ v )I * v = ( Nv -I * v ) m j=1 a j c j,1 ν j a j β v,j µ j + ν j µ j I * v µ j Nj + a j β v,j I * v -(µ v + δ v )I * v .
(2.9)

Since the equilibrium is endemic I * v > 0, Equation (2.9) is satisfied if, and only if:

0 = Nv -I * v µ v + δ v m j=1 c j,1 ν j a 2 j β v,j µ j + ν j µ j µ j Nj + a j β v,j I * v -1 := f (I * v ).
(2.10)

The function f (I * v ) is decreasing and f (0) = R 2 0 (m) -1. Moreover, f ( Nv ) = -1 < 0. Hence, the equation f (I * v ) has a unique positive root in (0, Nv ) if, and only if f (0) > 0, that is, whenever

R 2 0 (m) > 1.
Since the other components S * j , E * j , I * i,j , for j = 1, 2, . . . , m and i = 1, 2, . . . , n, are uniquely expressed in terms of I * v , we conclude that there is a unique endemic equilibrium in Int(Ω) whenever R 2 0 (m) > 1.

Theorem 2.3. The endemic equilibrium is globally asymptotically stable on IR m(n+2)+1 .

Proof.

We consider the following Lyapunov function

V = m j=1 A j0 S j -S * j log S j + E j -E * j log E j + m j=1 n i=1 A j,i I j,i -I * j,i log I j,i + V 1 (S v -S * v log S v ) + V 1 (I v -I * v log I v ) (2.11)
where the coefficients V 1 , A j0 , A ji , for i = 1, 2, . . . , n and j = 1, 2, . . . , m are positive to be determined later. The function V is thus positive-definite.

The derivative of V along trajectories is:

V = -V 1 (µ v + δ v )I v + A v + V 1 m j=1 a j S * v Nj n i=1 β j,i I * j,i 2 - S * v S v - S v S * v I j,i I * j,i I * v I v + a j S * v n i=1 β j,i I j,i Nj + m j=1 A hj + A j0 a j β v,j S * j I * v Nj 2 - S * j S j - S j S * j I v I * v E * j E j + A j0 a j β v,j S * j I v Nj -A j0 (µ j + ν j )E j + A j1 ν j E j -α j,1 I j,1 -A j1 ν j E j I * j,1 I j,1 + A j1 α j,1 I * j,1 + n i=2 A ji γ j,i-1 I j,i-1 -α j,i I j,i -A ji γ j,i-1 I j,i-1 I * j,i I j,i + A ji α j,i I * j,i , (2.12)
where

A v = (µ v + δ v )S * v 2 - S * v S v - S v S * v and A hj = µ h S * j 2 - S * j S j - S j S * j .
We choose the coefficients A ji in a similar fashion as in the one host case, that is,

A j0 (µ j + ν j )E * j = A j0 a j β v,j S * j I * v N j = V 1 a j S * v Nj n i=1 β j,i I * j,i for all j = 1, 2, . . . , m, (2.13)
and

A j1 = 1 ν j E * j V 1 a j S * v Nj n i=1 β j,i I * j,i
for all j = 1, 2, . . . , m.

Also,

A ji = n k=i β j,k I * j,k α j,i I * j,i V 1 a j β j,v S * v Nj for all 2 ≤ i ≤ n -1, 1 ≤ j ≤ m, and
α j,n A jn = V 1 a j β j,v S * v Nj for all 1 ≤ j ≤ m.
We obtain

V = V 1 m j=1 a j β j,1 S * v Nj I * j,1 5 - S * v S v - S v S * v I j,1 I * j,1 I * v I v - S * j S j - S j S * j I v I * v E * j E j - E j E * j I * j,1 I j,1 +V 1 m j=1 a j S * v Nj n i=2 β j,i I * j,i 4 + i - S * v S v - S v S * v I j,i I * j,i I * v I v - S * j S j - S j S * j I v I * v E * j E j - E j E * j I * j,1 I j,1 - i k=2 I j,k-1 I * j,k-1 I * j,k I j,k +A v + m j=1 A hj .
The function V is therefore negative-definite. Since V is positive-definite, the global stability of the interior endemic equilibrium on IR m(n+2)+1 follows from the Lyapunov stability theorem.

The result of global stability of the endemic equilibrium generalizes [START_REF] Cruz-Pacheco | Control measures for chagas disease[END_REF][START_REF]Multi-species interactions in west nile virus infection[END_REF] for an arbitrary number of hosts. Indeed, in [START_REF] Cruz-Pacheco | Control measures for chagas disease[END_REF], Cruz-Pacheco et al. considered and SIR-SI with two hosts in their model capturing Chagas' disease. The authors also considered the case where one host (humans) has two stages of infections, namely acute and chronic. Our model considers an arbitrary number of hosts and and host's infectious stages. Similarly, [START_REF]Multi-species interactions in west nile virus infection[END_REF] considered a model describing NWV dynamics with an arbitrary number of hosts, but the result of global stability has been done only for the two hosts case.

Remark 2.2.

We have proved that the multi-host model (2.3) satisfies the sharp-threshold property, that is, the disease dies out if R 2 0 (m) ≤ 1 and persists otherwise. However, R 2 0 (m) = m j=1 R 2 0,j , where R 2 0,j is the basic reproduction number in presence of Host j only. Hence,

R 2 0 (m) ≤ 1 =⇒ R 2 0,j ≤ 1 for all j = 1, 2, . . . , m.
Thus, the eradication of the the multi-host disease in one host is subject to the eradication of the disease from all hosts. However, the multi-host system could have a basic reproduction number great than one even though the reproduction number of each isolated host-vector is less than one (see Fig. 2 and Fig. 3).

In the next subsection, we provide a set of simulations that illustrate our analytical results and showcase that control strategies that target host-vector system in isolation may not be sufficient to eradicate VBZ.

Simulations

To illustrate the previous theoretical results, we run a set of simulations of a two-host system with four stages of infection, that is m = 2 and n = 4. Unless otherwise stated, the baseline parameters for the simulations are those given in Table 3. For these values, the overall basic reproduction number is R 2 0 (2) = 1.50. Figure 3 shows that the trajectories of the solutions representing the infected populations for all stages for both hosts tend to a nonzero endemic equilibrium as predicted by Theorems 2.2 and 2.3. However, if the two hosts are are considered separately, the system behave as two isolated host-vector systems, with the same parameter values, the corresponding basic reproduction numbers are R 2 0,1 = 0.66 and R 2 0,2 = 0.84 and therefore the disease dies out from both hosts and the vector population as portrayed in Fig. 2 and stated in Theorem 1.1. Figures 2 and3 exemplify also Remark 2.2, in which it is noted that controlling VBZ in each species in isolation, does not lead necessarily a complete control of VBZ. A coordinated effort in all host species and vector population is needed to steer the overall basic reproduction number below unity, and therefore eliminating the disease in all host species and vectors. Also, the simulations capture that host diversity increases the prevalence of the disease, leading to the "amplification effect." By choosing the values of βs lower than those given in Table 3 will lead to the overall basic reproduction number R 2 0 (2) ≤ 1 and the simulations show that all infected population tend to zero as proved in Theorem 2.1. Alternatively, if the parameters are chosen such that R 2 0,1 > 1 and R 2 0,2 > 1, the trajectories of the isolated hostvector systems converge to an endemic equilibria as expected (Theorem 1.2). We decided not to displays the figures of these cases.

Conclusion and discussions

Mathematical models of zoonoses have often been investigated when the transmission is direct, in which case, these models could be seen as multi-group models that well investigated in the literature. However, most zoonotic pathogens are spread through arthropod vectors such as mosquitoes, flies, fleas, ticks, etc., that transmit pathogens across species.

We considered a general host-vector SEIR-SI model where host's infectious state is subdivided into n classes, each of which has a different infectiosity to the vector. We derived the basic reproduction number R 0 . Our results show that the dynamics of the model is robust and is tied to R 0 . The disease-free equilibrium is globally asymptotically stable (GAS) if R 0 ≤ 1 and unstable otherwise (Theorem 1.1). A unique endemic equilibrium exists and is GAS whenever R 0 > 1 (Theorem 1.2).

We extended the host-vector staged progression model to m host species interacting with an arthropod vector to obtain a general framework in modeling a class of zoonoses. The basic reproduction number R 2 0 (m) for the system with m hosts is derived and happened to be the sum of basic reproduction number of the host-vector systems when a unique host is interacting with the arthropod vector, that is, R 2 0,j , for j = 1, 2, . . . , m. We proved that the system with m hosts remains robust as its asymptotic behavior is solely determined by R 0 (m). The disease dies out is above unity (Theorems 2.2 and 2.3).

However, the disease might persist in all hosts, that is, the global R 0 (m) could be greater that unity even if R 0,j < 1 for j = 1, 2 . . . , m, that is, for all hosts. Therefore, host heterogeneity favors the spread of the infection, and thus causes the so-called "amplification effect" [START_REF] Miller | The effects of host diversity on vector-borne disease: the conditions under which diversity will amplify or dilute the disease risk[END_REF] and opposed to the "dilution effect" [START_REF] Miller | The effects of host diversity on vector-borne disease: the conditions under which diversity will amplify or dilute the disease risk[END_REF][START_REF] Ostfeld | Effects of host diversity on infectious disease[END_REF][START_REF] Salkeld | A meta-analysis suggesting that the relationship between biodiversity and risk of zoonotic pathogen transmission is idiosyncratic[END_REF], for which the increase in host diversity could potentially decrease or drive out the pathogens [START_REF] Ostfeld | Effects of host diversity on infectious disease[END_REF]. As R 2 0 (m) = m j=1 R 2 0,j and considering the GAS of the unique endemic equilibrium, we deduced that the persistence of the disease in one host-vector system is sufficient to ensure persistence in all hosts. The overall mitigation or elimination of vector-borne zoonoses requires control strategies that bring the basic reproduction numbers in all hosts significantly below unity, and therefore reduces the disease burden in all hosts. Hence, this paper encompasses and reinforces the One Health [START_REF] Health | [END_REF] concept, for which human, animal and ecosystem health are ecologically interconnected [START_REF]One health, One Health Commission[END_REF].

Figure 1 :

 1 Figure 1: Flow diagram of Model 1.4.

  14) where b T = (b 1 , b 2 , . . . , b n ) and à is the submatrix A obtained by removing the first row and first column. We choose b and b 0 that cancel the coefficients of I h and E h , respectively. That is, b

Figure 2 :

 2 Figure 2: Simulation of two separate systems with one host each and four stages of infection. The R 20,j values are 0.66 for the host 1 population and 0.84 for the host 2 population. As expected, the infected populations go to zero.

Figure 3 :

 3 Figure3: Simulation of a two-host system where there are four stages of infection for each host. The R 2 0,j values are 0.66 for the host 1 population and 0.84 for the host 2 population, giving a combined system R 2 0 (2) = 1.50. The infected populations eventually go to the nonzero endemic equilibrium.

Table 1 :

 1 Description of the parameters used in System (1.4).

	Parameters	Description
	Λ h	Recruitment of the host
	Λ v	Recruitment of vectors
	a	Bitting rate
	µ h	Host's death rate
	β v,h	Host's infectiousness to mosquitoes per biting
	β i	Vector's infectiousness to host at stage i per biting
	ν h	Host's rate at which the exposed individuals become infected
	η i	Per capita recovery rate of an infected host at stage i
	γ i	Host's per capita progression rate from stage i to i + 1
	µ v	Vectors' natural mortality rate
	δ

v Vectors' control-induced mortality rate Remark 1.1.

  and c 1 and c 2 are defined by the relations (1.8).

		Hence, relations
	(1.11) and (1.12) imply that (S * h , I * h , I * v ) T	0 if and only if R 2 0 > 1. The next theorem gives the
	asymptotic behavior of this equilibrium whenever it exists.
	Theorem 1.2.

Table 2 :

 2 Description of the parameters used in System (2.2).

	Parameters	Description
	Λ j	Recruitment of Host j
	Λ v	Recruitment of vectors
	a j	Bitting rate on Host j
	µ j	Host j's death rate
	β v,j	Infectiousness of Host j to mosquitoes per biting
	β j,i	Vector's Infectiousness to host j at stage i per biting
	ν

j Host j's rate at which the exposed individuals become infected

η j,i

Per capita recovery rate for Host j at stage i γ j,i Host j's per capita progression rate from stage i to i + 1

Table 3 :

 3 Parameters used in the simulations.

	Parameter Host 1 Host 2
	Λ j	100	100
	Λ v	10000 10000
	β v,j	0.2	0.2
	µ j	1/28875 1/28875
	µ j	1/15	1/15
	ν j	1/15	1/15
	a j	0.37	0.37
	δ v	0.001	0.001
	β j,1	0.2	0.2
	β j,2	0.3	0.3
	β j,3	0.4	0.4
	β j,4	0.5	0.5
	η j,1	0.01	0.01
	η j,2	0.007	0.007
	η j,3	0.004	0.004
	η j,4	0.001	0.001
	γ j,1	0.125	0.125
	γ j,2	0.01	0.0875
	γ j,3	0.05	0.05
	γ j,4	0.025 0.0125
	R 2 0,j	0.66	0.84
	R 2 0 (2)	1.50

from all hosts if R 0 (m) ≤ 1 (Theorem 2.1) and persists in all hosts if R 0 (m) > 1. Moreover, we proved that a unique endemic equilibrium exists and is globally asymptotically stable if R 0 (m)