Nearly two-thirds of all known human infectious diseases (ID) are caused by zoonotic pathogens which are transmissible from one host species (humans and vertebrate animals) to another [25, 38] , and therefore multi-host. Moreover, 75% of emerging and re-emerging infectious diseases are classified as zoonoses and constitute a major public health problem around the world, responsible for over one million death and hundreds of millions of cases each year [42] .
Furthermore, it is estimated that zoonoses cause over 20 billion and 200 billion USD of direct and indirect economic burden across the world respectively [25, 42] . Although the morbidity and mortality of most ID decreased, the incidence of zoonoses have increased [26] . Therefore, understanding the dynamics of zoonoses by systematic modeling and analyzing in order to control and mitigate these scourges should be a worldwide priority.
Multi-host infectious diseases include Lyme disease, tick-borne relapsing fever (TBRF), West Nile virus (WNV), Chagas disease, type A influenzas, Rift Valley fever, severe acute respiratory syndrome (SARS), etc. However, 40% of multi-host ID are vector-borne [22] . That is, blood-sucking arthropod vectors such as ticks, mosquitoes, fleas and sandflies play the role of connecting multiple hosts while potentially infecting them and getting infected during this process. Thus, formulating the dynamics of zoonoses requires understanding the ecology of all involved host species and their interactions with vector species that carry the pathogen between hosts. These complexities have made a comprehensive study of zoonoses very challenging, making mathematical models of zoonoses scarce and mainly focused on WNV [8] or directly transmitted zoonoses [2, 3, 28] . Recent research has been focused on understanding the dynamics and control of vector-borne zoonoses with one vector and two hosts (see [32, 36] for WNV, [11] for Chagas' disease, and [23, 31] for TBRF). Also, authors in [6] proposed a class of vector-borne zoonoses with an arbitrary number of hosts and one vector and provided the complete global dynamics of equilibria.
However, another layer of complexity regarding the ecology of zoonoses consists of sometimes different arthropod vector species, or the same arthropod vectors but different genera, are responsible of transmiting the same pathogen to a number of different hosts. For instance, over 65 different mosquito species transmit WNV to a number of hosts including humans, mammals and many species of birds [17, 24] . Another example of a highly complex zoonosis is the eastern equine encephalitis virus (EEEV). Indeed, as illustrated in [10, 29] , the main vector of EEEV is the mosquito Culiseta melanura, which feeds exclusively on birds, and thus infects birds only. These birds infect, in turn, other mosquito species, which then bite humans, therefore creating a direct chain of transmission of the pathogen between a host and vector species with no direct link and/or transmissibility. Hence, as stated in [29] , the elimination or mitigation of zoonoses requires breaking the multiple transmission cycles corresponding to each potential host and vector species.
Therefore, a better understanding of zoonoses requires taking into account the ecology of all host species (including dead-end hosts -hosts that do not contribute to further transmission of the pathogen) and all vector species, along with the epidemiology of the zoonosis within the before-mentioned species. Moreover, developing general theories for the role of intermediate hosts in pathogen emergence is one of the nine challenges in modeling the emergence of novel pathogens according to [27] . The goal of this paper is twofold.
The first goal of this paper is to derive a class of models that describes the interaction between m host species and p vector species where the latter differ in their propensity to acquire or infect the pathogen from the former. Moreover, we structure each host's infectious stage into n different "ages" or classes where they infect each vector species at different rates. The proposed model provides a general framework in modeling zoonoses as it takes into account the complex patterns and multi-faceted host species dynamics along with their interactions with vector species. The derived class of models handles multiple levels of organizations including the case where some host or vector species have different epidemiological structures than others. Our modeling framework offers a plethora of possible scenarios and could be applied to study specific cases of zoonoses, and thus hopes to provide a forum that could help guide decisions on control efforts to mitigate and/or eradicate some zoonoses.
The second goal of the paper is to study the global asymptotic behavior of the proposed system. Particularly, we derive conditions under which the disease dies out or persists in all host species and vector species. A key threshold quantity happens to be the basic reproduction number R 2 0 (m, n, p) for the general system with m host, n stages, and p vector species. We prove that the disease-free equilibrium is globally asymptotically stable (GAS) whenever R 2 0 (m, n, p) < 1. The proof of the uniqueness of an endemic equilibrium (EE) for large epidemic systems is known to be challenging. Indeed, the uniqueness of an EE may not hold for multi-group directly transmitted diseases [35, 21] and sexually transmitted diseases [19] . For our system, we will prove indeed that it has a unique strongly EE (in the sense of Thieme [39] ) under the assumption that the host-vector network configuration is irreducible and R 2 0 (m, n, p) > 1. To do so, we transformed the system at equilibrium into an auxiliary dynamical system and showed that the equilibrium of the newly crafted system, is globally attractive under the pre-cited hypotheses. Moreover, we will prove that the strongly EE is GAS whenever it exists. The latter relies on a carefully constructed Lyapunov function and elements of graph theory. The authors in [14, 15] first used tools of graph theory to study the GAS of the EE for multi-group models for directed transmitted diseases. To do so, they derived a Lyapunov function for the multi-group system from that of a single epidemic system. Also, in [33, 34] , the authors used these tools to investigate for stage-progression models for water-borne diseases and the authors in [20] successfully generalized the approach to a multi-group vector-borne SIR − SI system. For our model, the arbitrary number of host species, stages, and vector species (all potentially different) present specificities that make the approach employed the after-mentioned papers impractical. For instance, the construction of the Lyapunov function from that of one group does not apply in our case.
The paper is outlined as follows: Section 1 is devoted to the derivation of the model. Section 2 lays out the basic the properties of the model and determines the associated basic reproduction number. Section 3 provides a complete analysis of the model, and Section 4 is devoted to concluding remarks and discussions.
We consider the dynamics of a disease transmitted by the interplay between m host species and p arthropod vectors. We assume that the disease follow an SEI n R − SI structure where n designates the number of infectious stages in the evolution of the disease within each host species. Let N i represents the total population of each host species i (1 ≤ i ≤ m) and let S i , E i , I l,i , and R i denote the susceptible, latent, infectious at stage l (1 ≤ l ≤ n), and recovered populations of host species i, respectively. Let N v,j be the total population of arthropod vectors of species j, (1 ≤ j ≤ p), each of which is composed by susceptible vectors S v,j and infectious vectors I v,j . The susceptible populations of host species i are generated through a constant recruitment Λ i , subjected to a natural death rate of µ i and could be infected after being bitten by an infectious vector of any species j. After being infected, the susceptible populations of Host i become latent, who then become infectious after an incubation period of 1/ν i . The infectiosity period of host species is structured into n stages, each characterizing the level of parasitemia in the corresponding host. At each stage l, the infectious of host i recover at a rate η l,i , progressing to the next stage of infection at a rate γ l,i or naturally die at a rate µ i . The susceptible arthropod vectors of species j, S v,j , could be infected by all infectious hosts, of any species and of any stage, at a different rates. Moreover, they are replenished at a constant rate Λ v,j and die either by natural death, at rate µ v,j , or due to control strategies specific for each vector of species j, at a rate δ v,j .
The overall multi-host, multi-stage and multi-vector model is captured by the following system:
A schematic description of the model is captured by Figure 1 and the parameters are described in Section 1. The total population of host species and vector species are asymptotically constant as their dynamics are given byṄ
Hence, by using the theory of asymptotic systems [9, 41] and by denoting the limits of N i and N v,j by
µ v,j +δ v,j , System (1) is asymptotically equivalent 
Vectors' natural mortality rates;
Vectors' control-induced mortality rates. to its limit system. Moreover, in this case, it could be written in a compact form as: . . , E m ] T , I k = [I k,1 , I k,2 , . . . , I k,m ] T is the vector of infectious at stage k (1 ≤ k ≤ n) for all hosts species. Also, α k,i = µ k,i + η k,i + γ k,i and γ n,i = 0 for all i. For vectors, . . , I v,p ] T denote the vectors of susceptible and infected, respectively. The matrices A, B ⋄ and B l are given by
A represents the biting/landing rates matrix, B ⋄ captures the infectiousness
Host 2 S2 E2 I1,2 I2,2 . . .
Host m of host species to vector species; and the matrices B k represents the infectiousness of vector species from infected host species at stage k. Model (2) describes the dynamics of a Multi-Host and Multi-Vector SEI n R − SI model where the infectious stage in each host is composed of n stages. The model is flexible and could be adapted to cases where different hosts could have different epidemiological structures with respect to the infection. This extends the work in [6, 23, 31] which consider the dynamics zoonoses and one vector species. Our model extend also the multi-host and multi-vector in [12] which explores and SIR-SI type of host-vector interaction. Particularly, our model consists of SEI n R where the infection of each host is stratified into n infectious classes and each of these classes infect susceptible vectors at different rates. This extension captures a more realistic aspect of infection in the interactions between hosts and vectors. Moreover, although we considered an SI structure for the vectors for simplicity, the results of this paper are also valid if the vector species's dynamics follows an SEI structure. Indeed, for some vectors such as mosquitoes, their incubation period is often nearly two weeks, which is on the same scale as their lifespan, making an SEI model more suited for their dynamics.
The set
where • denotes the Hadamard product, is a compact attracting positively invariant for System (2) . Therefore, the solutions of System (2) are biologically substantiated. The trivial equilibrium of System (2) is the disease-free equilibrium (D.F.E) and is given by
In the following lemma, we derive the basic reproduction number for Model (2) following [13, 40] .
Lemma 2. 1 . The basic reproduction number of Model (2) is given by:
where ρ(.) denotes the spectral radius operator.
Following the next generation method [13, 40] , the system composed of the infected variables in (2) could be decomposed as the sum of two columns vectors F (E, I) and V(E, I) where F represents new infections in each host and arthropod species and V that of transitions between classes. By letting F = DF (E, I)| E 0 and V = DV(E, I)| E 0 , the Jacobian of F and V, evaluated at the DFE, we obtain :
It is worth noticing that the matrix B has m × n columns and so, the last row of F has mn + m + p = m(n + 1) + p columns, in accordance with the first row of F , as A • B ⋄ has p columns. The matrix V is Metzler (positive off-diagonal entries) and therefore −V −1 ≥ 0 [4] . Given the particular shape of the matrix V , the matrix-wise entries of −V −1 are given by the induction relationship:
where α 0 = µ + ν and γ 0 = ν. Hence, we can deduce that
Hence, the next generation matrix is given by:
The basic reproduction number is the spectral of the next generation matrix. Hence, by denoting R 0 (m, n, p) the basic reproduction number of System (2) with n hosts, n infectious stages (in each host) and p vectors, we have:
Therefore,
This proves our claim.
In the next section, we study the asymptotic properties of System (2). Particularly, we show that the dynamics of the system is completely determined by R 2 0 (m, n, p) under certain conditions.
3 Global stability of equilibria 3 
The disease-free equilibrium (DFE) always exists and is in Ω. The following theorem gives conditions under which it is globally asymptotically stable (GAS).
Theorem 3. 1 . The DFE is GAS whenever R 2 0 (m, n, p) < 1.
The proof consists of proving that all infected variables of System (1) converge to zero and using the local stability of the DFE when R 2 0 (m, n, p) < 1 to conclude the GAS of the diseasefree steady state. Considering that diag(S h ) ≤ diag(N h ) anḋ
we obtain, from System (2), that
where F and V are the matrices generated in the next generation method. Since F is a nonnegative matrix and V is a Metzler matrix, we have (see [4] ),
is the stability modulus of F + V . Hence, the trajectories of the auxiliary system whose RHS is that of (3) converge to zero whenever R 0 (m, n, p) = ρ(−F V −1 ) < 1.
Since all the variables are positive, by the comparison theorem [37] , we conclude that lim t→∞ E = lim t→∞ I 1 = · · · = lim t→∞ I n = 0 m and lim
Moreover, this implies that:
where ⊘ denotes the Hadamard division. Relations (4) and (5) imply the attractivity of the DFE. Moreover, by [13, 40] , the DFE is locally asymptotically stable whenever R 2 0 (m, n, p) < 1. We conclude thus that the DFE is GAS whenever R 2 0 (m, n, p) < 1. Remark 3.1. In Theorem 3.1, it is worthwhile noting that no hypothesis on the irreducibility of the connectivity matrix is needed. This is important as this hypothesis is customarily used [5, 14, 20, 33] to obtain a Lyapunov function in the form of V = c T I where c is the left eigenvector of the next generation matrix.
In this subsection, we explore the asymptotic behavior of System (2) when R 2 0 (m, n, p) > 1. Particularly, we want to obtain conditions for which the disease persists in all stages, of all hosts and all vectors. That is, an equilibrium such that
Such interior equilibrium is also called strongly endemic [39] . The existence of such equilibrium is tied to the overall basic reproduction number and the connectivity between host and vector species. That is, the matrix
The following theorem gives the existence conditions of such an equilibrium.
Theorem 3.2. There exists a unique strongly endemic equilibrium for Model (2) whenever R 2 0 (m, n, p) > 1 and the Host-Vector connectivity configuration N is irreducible. Proof. An endemic equilibrium satisfies the relations:
In the following, we express all variables at equilibrium in terms of I * v . From (6), we express all I * l (1 ≤ l ≤ n) in terms of I * 1 , as follows:
From (7), it follows that I * l ≫ 0, for any 1 ≤ l ≤ n, if and only I * 1 ≫ 0. Moreover, we can also express I * 1 in terms of I * v and S * since
Using the equilibrium relation stemming from the equation of infected vectors (in 6), the relation S * v = N v − I * v , along with (7) and (8). we obtain
where
and α 0 = µ h + ν h and γ 0 = ν. Now, we want to express S in terms of I * v . That will make (9) in terms of I * v only. By using the equation of susceptible in (6), we obtain
Hence, Equation (9) leads to:
Thus, a unique strongly endemic equilibrium exists if and only if the vectorial equation (10) has a unique nonnegative solution. Moreover, notice that Equation (10) is satisfied if and only if I * v is an equilibrium of the auxiliary dynamical system:
Now, we will show that System (11) has unique equilibrium I * v ≫ 0 if the connectivity matrix N is irreducible and R 2 0 (m, n, p) > 1. To this end, we will use elements from monotone systems and particularly Hirsch's theorem [18] . Indeed, System (11) is monotone if and only if the vector field F (x) − diag(µ v + δ v )x is, or equivalently whenever the F ′ (x) is Metzler. That is, the off-diagonal elements of F ′ (x) are positive. We have:
since the first term is a diagonal matrix. It follows from (12) , that F ′ (x) is Metzler as long as
Hence, M and therefore F ′ (x) is Metzler. Moreover, since the Host-Vector configuration N is irreducible, F ′ (x) is irreducible. Thus, the auxiliary dynamical system (11) is strongly monotone. Moreover, it could seen that the map of F ′ (x) is monotonically decreasing. Since F (0 IR p ) − diag(µ v + δ v )0 IR p = 0 IR p and x ≤ N v , by Hirsch's theorem ( [18] , page 55), the trajectories of the auxiliary system (11) either tend to the origin, or else there is a unique equilibrium I * v ≫ 0 and all trajectories tend to I * v . The origin is unstable since, for G(
We have shown that if R 2 0 (m, n, p) > 1 and the host-vector configuration is irreducible, a unique solution I * v ≫ 0 of Equation (10) exists. Hence, using (6), (7), we deduce that E * ≫ 0 and I * l ≫ 0, for 1 ≤ l ≤ n. We conclude therefore that System (2) has a unique strongly endemic equilibrium whenever R 2 0 (m, n, p) > 1 and N is irreducible.
The following subsection consists of investigating the asymptotic properties of this unique strongly endemic equilibrium.
The strongly endemic equilibrium is GAS whenever it exists.
The coefficients v i are positive to be determined later. The function V is definite positive and the goal is show that its derivative along the trajectories of the multi-host, multi-vector system (2) is definite-negative. To ease the notations, let
Throughout this proof, we will be using the component-wise endemic relations, which could be written as:
The derivative of V along the trajectories of (2) is given by:
where
Given the expression of w ij , the terms in I v,j sum to zero. Similarly, using the expression of c 1,i , the terms in E i also sum to zero. Hence, what remains in (14) is:
Our first major claim is:
Indeed, the equality (⋆) claims that all linear terms in I l,i (1 ≤ l ≤ n) in (15) sum to zero.
By using the expressions of c k,i , we obtain: 
Showing (⋆⋆⋆) concludes the proof of (⋆). We notice that the relation (⋆⋆⋆) is satisfied if v i are the components of the solution of the linear systemBv = 0, with
where Moreover, we have: , which is exactly the right hand side of (⋆⋆⋆).
In summary, Relation (⋆), obtained through (⋆⋆) and (⋆⋆⋆) cancels linear terms in I l,i , for 1 ≤ l ≤ n, in the expression ofV given in relation (15) . Thus, the latter yields to:
By using endemic relations (13) and the expression of the coefficients c k,i , we have Substituting these expression and switching i by r in the last sum of (17), the expressionV becomes:
Also, notice that: Therefore, Equation (18) implies that:
Moreover, (19) implieṡ
Finally, by combining the second and third sums in (20) , we obtain:
The terms A v,j and A h,i are clearly definite-negative. Note that, taken separately, the two sums S 1 and S 2 in (21) are not definite negative. Hence, to finish the proof, we will prove that the sum S 1 + S 2 is definite-negative. To do so, we look at the coefficients of the sums from the graph-theoretical standpoint, following the same approach as [14, 15, 16, 20] . Indeed, let G(N ) be the directed graph that represents the connectivity N between the m hosts (including the l stages) and p vectors. Since the Host-Vector connectivity configuration N is irreducible, it follows that graph G(N ) is strongly connected. Recall that, for i = 1, 2, . . . , m, v i are components of the solution of the systemBv = 0, whereB is given by (16) . The matrix B is the so-called Laplacian matrix and its associated the graph G(B) is strongly connected sinceB is irreducible, since N is. Moreover, the solution ofBv = 0 is given by the Kirchhoff's matrix tree theorem [7, 30] as follows: where E(T ) is the set of all arcs in T . For our setup, an arc (i, l, i ′ , j) describes an infection arc that starts from Host i, at stage l directed to Host i ′ through Vector j. From a modeling standpoint, v i connects Host i to all vectors (throughβ i l,j ) and connect all vectors to all hosts but i (through w kj , with k = i). By using Cayley's formula [1, 30] , each v i = C ii is the sum of n m−1 p m−1 m m−2 terms, each of which is the product of m ′ − 1 w kjβ i l,j with k = i. These w kjβ i l,j represent the weight of each spanning tree T i , rooted at Host i.
, we investigate what each terms of v i w ijβ r l,j represents in terms of the graph G(B) and its spanning trees. Indeed, each term in v i w ijβ r l,j for all i, j, l, r is a weight of a unicyclic graph of a particular length, obtained by adding an arc (r, l, i, j) to a directed tree rooted at Host i, T ∈ T i . The obtained unicyclic graph Q has a unique cycle CQ of length 1 ≤ d ≤ m.
We group all isomorphic cycles as they will have the same coefficients. Hence, we conclude that S 1 + S 2 is the sum over all unicyclic graphs as S 1 + S 2 = Q S Q , where
Similarly, the product of the non constant terms in the last sum (22) is equal to 1 as CQ 3 is a cycle. Hence, by the arithmetic-geometric mean, S Q is definite-negative for each unicyclic graph Q. It follows from (21), thatV is definite-negative. Thus, the global stability of the EE follows from the Lyapunov stability theorem. This ends the proof.
Although the Lyapunov function obtained in Theorem 3.3 is somehow similar to those obtained in [14, 15] for multi-group models, they are structurally different, in the sense that it is not linear combination of Lyapunov functions of one-group (or one host, multiple vectors). Indeed, the orbital derivative of Lyapunov function V = m i=1 v i V i where is V i is the Lyapunov function for one host and multiple vectors, is not definite-definite negative for the coefficient v = (v 1 , v 2 , . . . , v m ) T .
In this paper, we have considered an SI model structure for the dynamics of the vectors. However, the results obtained here are valid for SEI type of structure. Indeed, in this case, the following Lyapunov function
, c 1,i = ν i + µ i ν i . This Lyapunov function has its derivative along the trajectories of the SEI n R − SEI system, definite negative. The corresponding coefficients v i are determined in the same fashion as the SEI n R − SI case.
In this paper, we formulated a multi-host, multi-stage and multi-vector epidemic model that describes the evolution of a class of zoonoses in which the pathogen is shared by multiple host species and the transmission occurs through the biting or landing of an arthropod vector. The proposed model improves those of [6, 23, 31] -by incorporating multiple arthropod vector species-and [11, 12] -by incorporating multiple hosts within multiple stages in each hosts' infectious class and the heterogeneous nature of the interactions between these host species and multiple arthropod vector species. We computed the basic reproduction number of the general system R 2 0 (m, n, p), for m hosts, n hosts' infectious stages and p vectors species. We proved that the disease free equilibrium is globally asymptotically stable whenever R 2 0 (m, n, p) < 1 (Theorem 3.1). Under the assumption the host-vector network configuration is irreducible, we proved that there exists a unique "strongly" endemic equilibrium as long as R 2 0 (m, n, p) > 1 and that, it is GAS whenever it exists. This result is new and improves previous results of [6] , for which the global result for multiple hosts and one vector is given, and of [12] in which a case of multi-species, multi-vector is considered but no stability results were given.
The global stability of the strongly equilibrium relies on a carefully constructed Lyapunov functions and tools of graph theory,à la [14, 15, 20] . The uniqueness and the global stability of the strongly endemic equilibrium requires the irreducibility of the host-vector network (Theorem 3.2 and Theorem 3.3), leading to the conclusion that the disease either dies out or persists in all hosts and all vectors. That is, controlling the disease requires an intervention in all hosts. This is close to impossible, given that some hosts are not even known. Furthermore, in some cases the natural habitats of some hosts and vectors are so apart that it is unlikely there is a direct transmission between these hosts and vector for an infection to take place. This could collapse the irreducibility of the hosts-vectors configuration. Thus, it is important to investigate the global asymptotic behavior of the solutions when the host-vector network configuration is reducible. This is the subject of a separate study and will be published elsewhere.
Other venues of expanding this work consist of considering different functional reproduction schemes for different vectors and/or hosts. Indeed, in this paper, although we have different reproduction (recruitment) rate for each hosts and vectors, they follow the same scheme, that is, they all consists of constant recruitments.