<figure><img src="image_1.png" /></figure>

HAL
open science

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

Derdei Bichara, Abderrahman Iggidr, Laura Smith

► To cite this version:

Derdei Bichara, Abderrahman Iggidr, Laura Smith. Multi-stage vector-borne zoonoses models: A global analysis. Bulletin of Mathematical Biology, 2018, 80 (7), pp.1810-1848. 10.1007/s11538-018-0435-1 . hal-01584474

HAL Id: hal-01584474

https://inria.hal.science/hal-01584474

Submitted on 8 Sep 2017

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.
Multi-stage vector-borne zoonoses models: A global analysis

Derdei Bichara$^{†}$, Abderrahman Iggidr$^{‡}$, Laura Smith$^{†}$

† Department of Mathematics & Center for Computational and Applied
Mathematics, California State University, Fullerton, CA 92831, USA

‡ Inria, Université de Lorraine, CNRS. Institut Elie Cartan de Lorraine,
UMR 7502. ISGMP Bat. A, Ile du Saulcy, 57045 Metz Cedex 01, France.

**Abstract**

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 \le 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.

**Mathematics Subject Classification:** 92D30, 34D23, 34D20, 34D40, 34A34.

**Keywords:** Vector-borne zoonoses, Stage progression, Multi-host, One Health, Amplification effect, Global stability, Nonlinear dynamical systems.

**Introduction**

Zoonoses, infectious diseases caused by pathogens transmissible under natural conditions from vertebrate animals to humans, account for 75% of emerging infectious diseases [46], with wildlife being an important source. Zoonotic pathogens represent 61% of all known human pathogens [46]. 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 [51] 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 major transmission routes [11]:
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 [26], 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 [32] 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 [1, 28] or along with a vector [39, 40, 41, 42].

While dynamical models capturing directly transmitted zoonoses, including food-borne and free-living pathogens, have received a substantial, yet insufficient, attention [2, 3, 6, 19, 23, 32] (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.) [7, 13, 14, 27, 44], 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 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 [27, 37] 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_0^2(m) = \sum_{j=1}^{m} R_{0,j}^2$, where $R_{0,j}$, for $j = 1, 2, \dots, 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) \le 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_0^2(m)$ in terms of $R_{0,j}^2$, for $j = 1, 2, \dots, 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 [12], which advocates for optimal health for
all species and the environment [35]. 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.”

# 1 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.

## 1.1 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, \dots, 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 + \sum_{i=1}^{n} 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 $\Lambda_h$ and reduced by a per-capita mortality rate of $\mu_h$ and through an infectious vector bite. The dynamics of the susceptible host is therefore given by:

$$ \dot{S}_h = \Lambda_h - b(N_h, N_v)\beta_{vh}S_h\frac{I_v}{N_v} - \mu_h S_h, \quad (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/\nu_h$. Hence, the dynamics of latent host is given by:

$$ \dot{E}_h = b(N_h, N_v)\beta_{vh}S_h\frac{I_v}{N_v} - (\mu_h + \nu_h)E_h. \quad (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 [29] and 5-15 days [10] 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 [17, 45]. For other tick-borne diseases, the incubation period varies from more than two weeks and up to nine weeks (see [9] 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 infectivity 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,

$$ \dot{I}_1 = \nu_h E_h - (\mu_h + \eta_1 + \gamma_1)I_1. $$

The infected host $I_1$ leave this stage either by natural death at a rate $\mu_h$, through recovery at a per capita rate $\eta_1$, or by progressing to the second stage of infectiousness at a per-capita rate of $\gamma_1$. Hence,

$$ \dot{I}_2 = \gamma_1 I_1 - (\mu_h + \eta_2 + \gamma_2)I_2. $$
Similarly, for $i = 3, 4, \dots, n-1$, the dynamics of the infected at stage $i$ is captured by:

$$ \dot{I}_i = \gamma_{i-1}I_{i-1} - (\mu_h + \eta_i + \gamma_i)I_i, $$

and

$$ \dot{I}_n = \gamma_{n-1}I_{n-1} - (\mu_h + \eta_n)I_n. $$

Hosts' stage-structure infectiousness in vector-borne zoonoses have long been recognized in literature [13, 24, 48]. 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 [45]. Moreover, according to [45], 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 [20, 22, 25], whereas the relapses is used more in ecology [27, 37, 45].

* 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 [18] after a seven year incubation period [50]. 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 [50].

* 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 $\eta_i$, the dynamics of the recovered are given by:

$$ \dot{R} = \sum_{i=1}^{n} \eta_i I_i - \mu_h R. $$

The susceptible vectors are recruited via a constant recruitment $\Lambda_v$ and die at natural death rate of $\mu_v$ or via control measures, such as Dichlorodiphenyltrichloroethane (DDT) (for mosquitoes, sandflies, etc.) and Acaricide (for ticks and mites) at a rate of $\delta_v$. The infection of susceptible vectors occurs after biting an infectious host. The total number of infectious hosts is $\sum_{i=1}^{n} 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 $\beta_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:
$$ \dot{S}_v = \Lambda_v - aS_v \sum_{i=1}^{n} \frac{\beta_i I_i}{N_h} - (\mu_v + \delta_v)S_v, $$

and

$$ \dot{I}_v = aS_v \sum_{i=1}^{n} \frac{\beta_i I_i}{N_h} - (\mu_v + \delta_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) = \frac{aN_v}{N_h}. \qquad (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:

$$ \left\{ \begin{aligned} \dot{S}_h &= \Lambda_h - a \beta_{vh} S_h \frac{I_v}{N_h} - \mu_h S_h \\ \dot{E}_h &= a \beta_{vh} S_h \frac{I_v}{N_h} - (\mu_h + \nu_h) E_h \\ \dot{I}_1 &= \nu_h E_h - (\mu_h + \eta_1 + \gamma_1) I_1 \\ \dot{I}_2 &= \gamma_1 I_1 - (\mu_h + \eta_2 + \gamma_2) I_2 \\ & \vdots \\ \dot{I}_{n-1} &= \gamma_{n-2} I_{n-2} - (\mu_h + \eta_{n-1} + \gamma_{n-1}) I_{n-1} \\ \dot{I}_n &= \gamma_{n-1} I_{n-1} - (\mu_h + \eta_n) I_n \\ \dot{R} &= \sum_{i=1}^{n} \eta_i I_i - \mu_h R \\ \dot{S}_v &= \Lambda_v - a S_v \sum_{i=1}^{n} \frac{\beta_i I_i}{N_h} - (\mu_v + \delta_v) S_v \\ \dot{I}_v &= a S_v \sum_{i=1}^{n} \frac{\beta_i I_i}{N_h} - (\mu_v + \delta_v) I_v \end{aligned} \right. \qquad (1.4) $$

The flow diagram of System (1.4) is represented in Fig. (1), 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 $\dot{N}_h = \Lambda_h - \mu_h N_h$ and $\dot{N}_v = \Lambda_v - (\mu_v + \delta_v) N_h$. Hence, we can deduce that:

$$ \lim_{t \to \infty} N_h = \frac{\Lambda_h}{\mu_h} := \bar{N}_h \quad \text{and} \quad \lim_{t \to \infty} N_v = \frac{\Lambda_v}{\mu_v + \delta_v} := \bar{N}_v. $$

By using the theory of asymptotically autonomous systems for triangular systems [8, 49], we can
<figure><img src="image_2.png" /><figcaption>Figure 1: Flow diagram of Model 1.4.</figcaption></figure>

see that System (1.4) is asymptotically equivalent to:

$$ 
\left\{
\begin{aligned}
\dot{S}_h &= \Lambda_h - a \beta_{vh} S_h \frac{I_v}{\bar{N}_h} - \mu_h S_h \\
\dot{E}_h &= a \beta_{vh} S_h \frac{I_v}{\bar{N}_h} - (\mu_h + \nu_h) E_h \\
\dot{I}_1 &= \nu_h E_h - (\mu_h + \eta_1 + \gamma_1) I_1 \\
\dot{I}_2 &= \gamma_1 I_1 - (\mu_h + \eta_2 + \gamma_2) I_2 \\
&\vdots \\
\dot{I}_{n-1} &= \gamma_{n-2} I_{n-2} - (\mu_h + \eta_{n-1} + \gamma_{n-1}) I_{n-1} \\
\dot{I}_n &= \gamma_{n-1} I_{n-1} - (\mu_h + \eta_n) I_n \\
\dot{I}_v &= a(\bar{N}_v - I_v) \sum_{i=1}^{n} \frac{\beta_i I_i}{\bar{N}_h} - (\mu_v + \delta_v) I_v
\end{aligned}
\right.
\quad (1.5)
$$

Variations of this model have been considered in the literature. Indeed in modeling tick-borne relapsing fever (TBRF), Johnson et al. in [27] and Palmer et al. [37] 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 [21, 22] 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 $\frac{1}{\nu} \to 0$, that is, when the mean latency period is 0, Model (1.5) is the limit of an SIR-SI model proposed in [27, 37]. Also, Park [38] 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.
<table><caption>Table 1: Description of the parameters used in System (1.4).</caption><thead><tr><th>Parameters</th><th>Description</th></tr></thead><tbody><tr><td>Λ<sub>h</sub></td><td>Recruitment of the host</td></tr><tr><td>Λ<sub>v</sub></td><td>Recruitment of vectors</td></tr><tr><td>a</td><td>Bitting rate</td></tr><tr><td>μ<sub>h</sub></td><td>Host's death rate</td></tr><tr><td>β<sub>v,h</sub></td><td>Host's infectiousness to mosquitoes per biting</td></tr><tr><td>β<sub>i</sub></td><td>Vector's infectiousness to host at stage i per biting</td></tr><tr><td>ν<sub>h</sub></td><td>Host's rate at which the exposed individuals become infected</td></tr><tr><td>η<sub>i</sub></td><td>Per capita recovery rate of an infected host at stage i</td></tr><tr><td>γ<sub>i</sub></td><td>Host's per capita progression rate from stage i to i + 1</td></tr><tr><td>μ<sub>v</sub></td><td>Vectors' natural mortality rate</td></tr><tr><td>δ<sub>v</sub></td><td>Vectors' control-induced mortality rate</td></tr></tbody></table>

**Remark 1.1.**

In [27], the authors assume different natural death rates in each of the infected classes but assume that $\mu_i = \frac{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 $\Lambda_h = \mu_S S_h + \mu_E E_h + \sum_{i=1}^n \mu_{h,i} I_i + \mu_{RR}$, 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:

$$ egin{cases} 
\dot{S}_h = \Lambda_h e_1 - a \beta_{vh} S_h \frac{I_v}{\bar{N}_h} - \mu_h S_h \\ 
\dot{I}_h = a \beta_{vh} S_h \frac{I_v}{\bar{N}_h} e_1 + A \mathbf{I}_h \\ 
\dot{I}_v = a \frac{\bar{N}_v - I_v}{N_h} \langle \begin{pmatrix} 0 \\ \beta \end{pmatrix} | \mathbf{I}_h \rangle - (\mu_v + \delta_v) I_v 
\end{cases} 
\quad (1.6) $$

where $\mathbf{I}_h = (E_h, I_1, \dots, I_n)^T$, $e_1 = (1, 0, \dots, 0)^T \in \mathbb{R}^{n+1}$, $\beta = (\beta_1, \beta_2, \dots, \beta_n)^T$, $\langle x | y \rangle$ is the canonical scalar product (here in $\mathbb{R}^{n+1}$), and $A$ is given by:

$$ A = \begin{pmatrix} -(\nu + \mu_h) & 0 & \dots & \dots & 0 \\ \nu & -\alpha_1 & \dots & \dots & 0 \\ 0 & \gamma_1 & -\alpha_2 & \dots & 0 \\ \vdots & \ddots & \ddots & \ddots & \vdots \\ 0 & \dots & \ddots & \gamma_{n-1} & -\alpha_n \end{pmatrix}, $$

where $\alpha_i = \gamma_i + \eta_i + \mu_h$, for $i = 1, 2, \dots, n-1$ and $\alpha_n = \eta_n + \mu_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

$$ \Omega = \left\{ (S_h, E_h, I_i, I_v) \in \mathbb{R}_+^{n+3} \mid S_h + E_h + \sum_{i=1}^n I_i \le \bar{N}_h, I_v \le \bar{N}_v \right\} $$

is a compact attracting positively invariant set for system (1.6).

The disease-free equilibrium is $E_0 = (\bar{S}_h, 0_{\mathbb{R}^{n+2}})$, where $\bar{S}_h = \frac{\Lambda}{\mu_h}$. This equilibrium always exists and it belongs to $\Omega$.

The infected compartments of System (1.6) are the sum of $\mathcal{F}$ and $\mathcal{V}$ where

$$ \mathcal{F}(\mathbf{I}_h, I_v) = \left( \begin{array}{c} a \beta_{vh} S_h \frac{I_v}{\bar{N}_h} e_1 \\ a \frac{\bar{N}_v - I_v}{N_h} \langle \begin{pmatrix} 0 \\ \boldsymbol{\beta} \end{pmatrix} \mid \mathbf{I}_h \rangle \end{array} \right) \quad \text{and} \quad \mathcal{V}(\mathbf{I}_h, I_v) = \left( \begin{array}{c} A \mathbf{I}_h \\ -(\mu_v + \delta_v) I_v \end{array} \right). $$

Let $F = D\mathcal{F}(\mathbf{I}_h, I_v)$ and $V = D\mathcal{V}(\mathbf{I}_h, I_v)$ be function-valued matrices evaluated at the DFE. We obtain:

$$ V = \begin{pmatrix} A & 0 \\ \mathbf{0}_{1,n+1} & -(\mu_v + \delta_v) \end{pmatrix} \quad \text{and} \quad F = \begin{pmatrix} \mathbf{0}_{n+1,n+1} & a\beta_{vh}e_1 \\ a\beta_{hv}\frac{\bar{N}_v}{\bar{N}_h}(0,\boldsymbol{\beta})^T & 0 \end{pmatrix}. $$

Given the particular form of $A$, it is not difficult to verify that:

$$ -A^{-1} = \begin{pmatrix} \frac{1}{\nu + \mu_h} & 0 & \dots & \dots & 0 \\ \frac{\nu}{\alpha_1(\nu + \mu_h)} & \frac{1}{\alpha_1} & \dots & \dots & 0 \\ \frac{\gamma_1 \nu}{\alpha_1 \alpha_2(\nu + \mu_h)} & \frac{\gamma_1}{\alpha_1 \alpha_2} & \frac{1}{\alpha_2} & \dots & 0 \\ \vdots & \ddots & \ddots & \ddots & \vdots \\ \frac{\gamma_1 \dots \gamma_{n-1} \nu}{\alpha_1 \dots \alpha_n(\nu + \mu_h)} & \dots & \dots & \frac{\gamma_{n-1}}{\alpha_{n-1} \alpha_n} & \frac{1}{\alpha_n} \end{pmatrix} \quad (1.7) $$

Hence, the next generation matrix is given by:

$$ -FV^{-1} = \begin{pmatrix} \mathbf{0}_{n+1,n+1} & \frac{a\beta_{vh}}{\mu_v + \delta} e_1 \\ -a \frac{\bar{N}_v}{\bar{N}_h} (0, \boldsymbol{\beta})^T A^{-1} & 0 \end{pmatrix}. $$

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 [15, 47]. Denoted by $\mathcal{R}_0$, it is the spectral radius of the next generation matrix $-FV^{-1}$. Therefore, for Model (1.6), the basic reproduction number is,

$$ \begin{aligned} \mathcal{R}_0^2 &= \frac{a^2 \beta_{vh} \bar{N}_v}{\mu_v + \delta} (0, \boldsymbol{\beta})^T (-A^{-1}) e_1 \\ &:= \frac{a^2 \beta_{vh} \bar{N}_v \nu}{\bar{N}_h (\mu_v + \delta) (\nu + \mu_h)} \sum_{i=1}^n \frac{\gamma_1 \gamma_2 \dots \gamma_{i-1}}{\alpha_1 \alpha_2 \dots \alpha_i} \beta_i \quad \text{with } \gamma_0 = 1. \end{aligned} $$

This quantity is intrinsically tied to the dynamics of the disease as we will discuss in the next subsection.
## 1.2 Global stability of the DFE

By the derivation of the basic reproduction number, the DFE is locally asymptotically stable if $\mathcal{R}_0 < 1$ and unstable if $\mathcal{R}_0 > 1$ [15, 47]. The following theorem provides a global result.

**Theorem 1.1.**

The DFE is globally asymptotically stable if $\mathcal{R}_0^2 \le 1$.

*Proof.*

The instability of the DFE when $\mathcal{R}_0 > 1$ follows from [47]. As for the global stability of the DFE when $\mathcal{R}_0 \le 1$, let us consider the candidate Lyapunov function:

$$ \mathcal{V}(\mathbf{I}_h, I_v) = \mathcal{L}(\mathbf{I}_h) + c_v I_v $$

where

$$ \mathcal{L}(\mathbf{I}_h) = c_0 E_h + \sum_{i=1}^{n} c_i I_i, $$

with $c_v = \frac{\bar{N}_h}{a\bar{N}_v}$ and the coefficients $c_i$ for $i = 1, 2, \dots, n-1$ are given by the following induction relationship:

$$ \left\{ \begin{array}{l} c_0 = \frac{\nu}{\nu + \mu_h} c_1 \\ c_1 = \frac{\beta_1}{\alpha_1} + \frac{\gamma_1}{\alpha_1 \alpha_2} \beta_2 + \dots + \frac{\gamma_1 \dots \gamma_{n-1}}{\alpha_1 \alpha_2 \dots \alpha_n} \beta_n, \\ c_{i+1} = \frac{1}{\gamma_i} (\alpha_i c_i - \beta_i), \quad \text{for } i = 1, 2, \dots, n-2, \\ c_n = \frac{\beta_n}{\alpha_n}. \end{array} \right. \qquad (1.8) $$

The derivative along the trajectories of the system is:

$$ \begin{aligned} \dot{\mathcal{V}}(\mathbf{I}_h, I_v) &= \dot{\mathcal{L}}(\mathbf{I}_h) + c_v \dot{I}_v \\ &= c_0 \dot{E}_h + \sum_{i=1}^{n} c_i \dot{I}_i + c_v \dot{I}_v \\ &= c_0 \left( a \beta_{vh} S_h \frac{I_v}{\bar{N}_h} - (\mu_h + \nu_h) E_h \right) + c_1 (\nu_h E_h - \alpha_1 I_1) + \sum_{i=2}^{n} c_i \dot{I}_i + c_v \dot{I}_v \\ &= c_0 a \beta_{vh} S_h \frac{I_v}{\bar{N}_h} - c_1 \alpha_1 I_1 + c_2 (\gamma_1 I_1 - \alpha_2 I_2) + \sum_{i=3}^{n} c_i \dot{I}_i + c_v \dot{I}_v \\ &= c_0 a \beta_{vh} S_h \frac{I_v}{\bar{N}_h} - c_1 \alpha_1 I_1 + (c_1 \alpha_1 - \beta_1) I_1 - c_2 \alpha_2 I_2 + \sum_{i=3}^{n} c_i \dot{I}_i + c_v \dot{I}_v, \end{aligned} $$

since $c_2\gamma_1 = \alpha_1c_1 - 1$. Therefore,

$$ \dot{\mathcal{V}}(\mathbf{I}_h, I_v) = c_0 a \beta_{vh} S_h \frac{I_v}{\bar{N}_h} - \beta_1 I_1 - c_2 \alpha_2 I_2 + \sum_{i=3}^{n} c_i \dot{I}_i + c_v \dot{I}_v. $$

Similarly, by using successively the induction relationship (1.8), we obtain,

$$ \dot{\mathcal{V}}(\mathbf{I}_h, I_v) = c_0 a \beta_{vh} S_h \frac{I_v}{\bar{N}_h} - \sum_{i=1}^{n} \beta_i I_i + c_v \left( a(N_v - I_v) \sum_{i=1}^{n} \frac{\beta_i I_i}{\bar{N}_h} - (\mu_v + \delta_v) I_v \right). \qquad (1.9) $$
Considering the expression of $c_v = \frac{\bar{N}_h}{a\bar{N}_v}$ and the inequality $S_h \le \bar{N}_h$, Equality (1.9) leads to:

$$
\begin{align*}
\dot{V}(\mathbf{I}_h, I_v) &= c_0 a \beta_{vh} \frac{S_h}{\bar{N}_h} I_v - I_v \sum_{i=1}^{n} \beta_i I_i - \frac{\bar{N}_h(\mu_v + \delta_v)}{a\bar{N}_v} I_v \\
&= \left( c_0 a \beta_{vh} \frac{S_h}{\bar{N}_h} - \frac{\bar{N}_h(\mu_v + \delta_v)}{a\bar{N}_v} \right) I_v - I_v \sum_{i=1}^{n} \beta_i I_i \\
&= \frac{\bar{N}_h(\mu_v + \delta_v)}{a\bar{N}_v} \left( \frac{a^2 \beta_{vh} \bar{N}_v c_0}{(\mu_v + \delta_v)\bar{N}_h} \frac{S_h}{\bar{N}_h} - 1 \right) I_v - I_v \sum_{i=1}^{n} \beta_i I_i \\
&= \frac{\bar{N}_h(\mu_v + \delta_v)}{a\bar{N}_v} \left( \mathcal{R}_0^2 \frac{S_h}{\bar{N}_h} - 1 \right) I_v - I_v \sum_{i=1}^{n} \beta_i I_i,
\end{align*}
$$

since $c_0 = \frac{\nu}{\nu + \mu_h} \sum_{i=1}^{n} \frac{\gamma_1 \gamma_2 \cdots \gamma_{i-1}}{\alpha_1 \alpha_2 \cdots \alpha_i} \beta_i$. Thus, we conclude that

$$ \dot{V}(\mathbf{I}_h, I_v) \le 0, $$

whenever $\mathcal{R}_0^2 \le 1$ as $S_h \le \bar{N}_h$. Moreover, $\dot{V}(\mathbf{I}_h, I_v) = 0$ if $I_v = 0$ or if $\mathcal{R}_0^2 = 1$, $S_h = \bar{N}_h$ and $\sum_{i=1}^{n} \beta_i I_i = 0$. Hence, we deduce that the largest set contained in $\{(\mathbf{I}_h, I_v) \in \Omega \mid \dot{V}(\mathbf{I}_h, I_v) = 0\}$ is reduced to the DFE. Since the set $\Omega$ is compact and positively invariant, by LaSalle's invariance principle [5, 31], the DFE is globally asymptotically stable on $\Omega$. The attractiveness of $\Omega$ makes the DFE globally asymptotically stable on $\mathbb{R}^{n+3}$. $\square$

## 1.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^*, \mathbf{I}_h^*, I_v^*)^T \gg 0$ be an endemic equilibrium for Model (1.6). Therefore $(S_h^*, \mathbf{I}_h^*, I_v^*)^T$ satisfies the endemic relations:

$$
\left\{
\begin{array}{l}
0 = \Lambda_h - a \beta_{vh} S_h^* \frac{I_v^*}{\bar{N}_h} - \mu_h S_h^* \\
0 = a \beta_{vh} S_h \frac{I_v^*}{\bar{N}_h} - (\mu_h + \nu_h) E_h^* \\
E_h^* = \frac{\alpha_1}{\nu_h} I_1^* \\
I_2^* = \frac{\gamma_1}{\alpha_2} I_1^* \\
I_3^* = \frac{\gamma_2}{\alpha_3} I_2^* = \frac{\gamma_1 \gamma_2}{\alpha_2 \alpha_3} I_1^* \\
\vdots \\
I_{n-1}^* = \frac{\gamma_{n-2}}{\alpha_{n-1}} I_{n-2}^* = \frac{\gamma_1 \gamma_2 \cdots \gamma_{n-2}}{\alpha_2 \alpha_3 \cdots \alpha_{n-1}} I_1^* \\
I_n^* = \frac{\gamma_{n-1}}{\alpha_n} I_{n-1}^* = \frac{\gamma_1 \gamma_2 \cdots \gamma_{n-1}}{\alpha_2 \alpha_3 \cdots \alpha_n} I_1^* \\
0 = a(\bar{N}_v - I_v^*) \sum_{i=1}^n \frac{\beta_i I_i^*}{\bar{N}_h} - (\mu_v + \delta_v) I_v^*
\end{array}
\right.
\quad (1.10)
$$

where $\alpha_i = \mu_i + \eta_i + \gamma_i$. It follows from (1.10) that there is an interior endemic equilibrium if $I_1^* > 0$. Moreover, we can deduce that:

$$
\begin{align}
\sum_{i=1}^{n} \beta_i I_i^* &= I_1^* \left( \beta_1 + \frac{\gamma_1}{\alpha_2} \beta_2 + \cdots + \frac{\gamma_1 \gamma_2 \cdots \gamma_{n-1}}{\alpha_2 \cdots \alpha_n} \beta_n \right) \nonumber \\
&= \alpha_1 c_1 I_1^*. \tag{1.11}
\end{align}
$$
The first equation of (1.10) implies that:

$$S_h^* = \frac{\mu_h \bar{N}_h}{\mu_h + a\beta_{vh} \frac{I_v^*}{\bar{N}_h}} \quad \text{and} \quad E_h^* = \frac{a\beta_{vh}}{\mu_h + \nu_h} \frac{\mu_h I_v^*}{\mu_h + a\beta_{vh} \frac{I_v^*}{\bar{N}_h}}.$$

We have

$$I_1^* = \frac{\mu_h(\mu_v + \delta_v)(\mathcal{R}_0^2 - 1)}{a\mu_h \bar{N}_h \xi + a^2 \beta_{vh} \bar{N}_v \xi}, \quad (1.12)$$

where $\xi = \beta_1 + \gamma_1 c_2 = \alpha_1 c_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 \gg 0$ if and only if $\mathcal{R}_0^2 > 1$. The next theorem gives the asymptotic behavior of this equilibrium whenever it exists.

**Theorem 1.2.**

The endemic equilibrium for the multi-host Model (1.6) is globally asymptotically stable whenever $\mathcal{R}_0^2 > 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) + \sum_{i=1}^{n} 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^*) + \sum_{i=1}^{n} b_i(I_i^* - I_i^* \log I_i^*)$$

and $b_i$, for $i = 0, 1, \dots, n$, and $v_1$ are positive constants. The function $W$ is therefore positive-definite 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

$$\bar{\beta}_i = a\beta_i \frac{1}{\bar{N}_h}, \quad \bar{\beta}_h = (\bar{\beta}_1, \dots, \bar{\beta}_n)^T \quad \text{and} \quad \bar{\beta}_v = a\beta_{vh} \frac{1}{\bar{N}_h}.$$

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

$$\begin{aligned} \dot{W} = b_0 \left(1 - \frac{S_h^*}{S_h}\right) (\Lambda_h - \bar{\beta}_v S_h I_v - \mu_h S_h) &+ b_0 \left(1 - \frac{E_h^*}{E_h}\right) (\bar{\beta}_v S_h I_v - (\mu_h + \nu_h) E_h) \\ &+ b_1 \left(1 - \frac{I_1^*}{I_1}\right) (\nu_h E_h - (\mu_h + \eta_1 + \gamma_1) I_1) + \sum_{i=2}^n b_i \left(1 - \frac{I_i^*}{I_i}\right) (\gamma_{i-1} I_{i-1} - (\mu_h + \eta_i + \gamma_i) I_i) \\ &+ v_1 \left(1 - \frac{S_v^*}{S_v}\right) (\Lambda_v - S_v \sum_{i=1}^n \bar{\beta}_i I_i - (\mu_v + \delta_v) S_v) \\ &+ v_1 \left(1 - \frac{I_v^*}{I_v}\right) (S_v \sum_{i=1}^n \bar{\beta}_i I_i - (\mu_v + \delta_v) I_v). \end{aligned} \quad (1.13)$$

By grouping all linear terms of (1.13) in $E_h$, $I_i$, for $i = 1, 2, \dots, n$ and $I_v$, and denoting it by $L$, we obtain after lengthy computations and arrangements:

$$L = (b^T \tilde{A} + v_1 S_v^* \bar{\beta}_h^T) \mathbf{I}_h + (b_1 \nu_h - b_0(\mu_h + \nu_h)) E_h + (b_0 \bar{\beta}_v S_h^* - v_1(\mu_v + \delta_v)) I_v \quad (1.14)$$

where $b^T = (b_1, b_2, \dots, b_n)$ and $\tilde{A}$ is the submatrix $A$ obtained by removing the first row and first column. We choose $b$ and $b_0$ that cancel the coefficients of $\mathbf{I}_h$ and $E_h$, respectively. That is,

$$\tilde{b} = -v_1 S_v^* \tilde{A}^{-T} \bar{\beta}_h \quad \text{and} \quad b_0 = \frac{\nu_h}{\mu_h + \nu_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:

$$
\begin{align*}
b_1 &= v_1 S_v^* \sum_{i=1}^{n} \frac{\gamma_1 \dots \gamma_{i-1} \bar{\beta}_i}{\alpha_1 \dots \alpha_i} \\
&= v_1 S_v^* \left( \frac{1}{\alpha_1} \bar{\beta}_1 + \frac{\gamma_1}{\alpha_1 \alpha_2} \bar{\beta}_2 + \dots + \frac{\gamma_1 \dots \gamma_{n-1}}{\alpha_1 \alpha_2 \dots \alpha_n} \bar{\beta}_n \right) \\
&= v_1 S_v^* \frac{1}{\alpha_1} \left( \bar{\beta}_1 + \frac{\gamma_1}{\alpha_2} \bar{\beta}_2 + \dots + \frac{\gamma_1 \dots \gamma_{n-1}}{\alpha_2 \dots \alpha_n} \bar{\beta}_n \right) \\
&= v_1 S_v^* \frac{1}{\alpha_1} \frac{\sum_{i=1}^{n} \bar{\beta}_i I_i^*}{I_1^*}. \tag{1.15}
\end{align*}
$$

Also, $b_0$ could clearly be expressed as

$$
\begin{align}
b_0 &= \frac{\nu_h}{\mu_h + \nu_h} b_1 \nonumber \\
    &= \frac{\nu_h}{\mu_h + \nu_h} v_1 S_v^* \frac{1}{\alpha_1} \frac{\sum_{i=1}^n \bar{\beta}_i I_i^*}{I_1^*}. \tag{1.16}
\end{align}
$$

The use of (1.16) cancels the coefficient of $I_v$ in (1.14). Indeed,

$$
\begin{align*}
b_0 \bar{\beta}_v S_h^* - v_1(\mu_v + \delta_v) &= \frac{\nu_h}{\mu_h + \nu_h} v_1 S_v^* \frac{1}{\alpha_1} \frac{\sum_{i=1}^n \bar{\beta}_i I_i^*}{I_1^*} \bar{\beta}_v S_h^* - v_1(\mu_v + \delta_v) \\
&= v_1 \frac{\nu_h}{(\mu_h + \nu_h) \alpha_1 I_1^*} S_v^* \sum_{i=1}^n \bar{\beta}_i I_i^* \bar{\beta}_v S_h^* - v_1(\mu_v + \delta_v) \\
&= v_1 \frac{\nu_h}{(\mu_h + \nu_h) \alpha_1 I_1^*} (\mu_v + \delta_v) I_v^* \bar{\beta}_v S_h^* - v_1(\mu_v + \delta_v) \\
&= v_1 \frac{\nu_h}{(\mu_h + \nu_h) \alpha_1 I_1^*} (\mu_v + \delta_v) (\mu_h + \nu_h) E_h^* - v_1(\mu_v + \delta_v) \\
&= v_1 \frac{\nu_h}{(\mu_h + \nu_h) \alpha_1 I_1^*} (\mu_v + \delta_v) (\mu_h + \nu_h) \frac{\alpha_1 I_1^*}{\nu_h} - v_1(\mu_v + \delta_v) \\
&= 0.
\end{align*}
$$

Therefore, with this choice for b and $b_0$, all linear terms in (1.13) cancel. The remaining terms in $\dot{W}$ are:

$$
\begin{align*}
\dot{W} &= v_1 \left(1 - \frac{S_v^*}{S_v}\right) \left(\Lambda_v - (\mu_v + \delta_v)S_v\right) + v_1(\mu_v + \delta_v)I_v^* - v_1 \frac{I_v^*}{I_v} S_v \sum_{i=1}^n \bar{\beta}_i I_i \\
&\quad + b_0 \left(1 - \frac{S_h^*}{S_h}\right) \left(\Lambda_h - \mu_h S_h\right) + b_0 \left( (\mu_h + \nu_h)E_h^* - \bar{\beta}_v S_h I_v \frac{E_h^*}{E_h} \right) + b_1(\mu_h + \eta_1 + \gamma_1)I_1^* \\
&\quad - b_1 \nu_h \frac{I_1^* E_h}{I_1} + \sum_{i=2}^{n-1} \left(b_i(\mu_h + \eta_i + \gamma_i)I_i^* - b_i \gamma_{i-1} \frac{I_i^* I_{i-1}}{I_i}\right) + b_n(\mu_h + \eta_n)I_n^* - b_n \gamma_{n-1} \frac{I_n^* I_{n-1}}{I_n} \\
&= v_1 (\mu_v + \delta_v) S_v^* \underbrace{\left(2 - \frac{S_v^*}{S_v} - \frac{S_v}{S_v^*}\right)}_{\mathcal{A}_v} + 2v_1(\mu_v + \delta_v)I_v^* - v_1(\mu_v + \delta_v) \frac{I_v^* S_v^*}{S_v} - v_1 \frac{I_v^*}{I_v} S_v \sum_{i=1}^n \bar{\beta}_i I_i \\
&\quad + b_0 \mu_h S_h^* \underbrace{\left(2 - \frac{S_h^*}{S_h} - \frac{S_h}{S_h^*}\right)}_{\mathcal{A}_h} + 2b_0(\mu_h + \nu_h)E_h^* - b_0(\mu_h + \nu_h) \frac{E_h^* S_h^*}{S_h} + b_0 \left(-\bar{\beta}_v S_h I_v \frac{E_h^*}{E_h}\right) \\
&\quad + b_1(\mu_h + \eta_1 + \gamma_1)I_1^* - b_1 \nu_h \frac{I_1^* E_h}{I_1} + \sum_{i=2}^{n-1} \left(b_i(\mu_h + \eta_i + \gamma_i)I_i^* - b_i \gamma_{i-1} \frac{I_i^* I_{i-1}}{I_i}\right) \\
&\quad + b_n(\mu_h + \eta_n)I_n^* - b_n \gamma_{n-1} \frac{I_n^* I_{n-1}}{I_n}.
\end{align*}
$$
By using the endemic relations and (1.15-1.16), $b_0$ and $b_1$ could be written as:

$$ 
\begin{aligned}
b_0 &= \frac{\nu_h}{\mu_h + \nu_h} b_1 \\
&= \frac{\nu_h}{\mu_h + \nu_h} v_1 S_v^* \frac{1}{\alpha_1} \frac{\sum_{i=1}^n \bar{\beta}_i I_i^*}{I_1^*} \\
&= \frac{v_1}{\mu_h + \nu_h} S_v^* \frac{\sum_{i=1}^n \bar{\beta}_i I_i^*}{E_h^*} \\
&= \frac{v_1}{\bar{\beta}_v} S_v^* \frac{\sum_{i=1}^n \bar{\beta}_i I_i^*}{S_h^* I_v^*}.
\end{aligned}
$$

Hence,

$$ 
\begin{aligned}
\dot{W} ={}& v_1 \mathcal{A}_v + 2v_1(\mu_v + \delta_v)I_v^* - v_1(\mu_v + \delta_v)\frac{I_v^* S_v^*}{S_v} - v_1 \frac{I_v^*}{I_v} S_v \sum_{i=1}^n \bar{\beta}_i I_i \\
& + b_0 \mathcal{A}_h + 2b_0(\mu_h + \nu_h)E_h^* - b_0(\mu_h + \nu_h)\frac{E_h^* S_h^*}{S_h} - b_0 \bar{\beta}_v S_h \frac{E_h^* I_v}{E_h} \\
& + b_1(\mu_h + \eta_1 + \gamma_1)I_1^* - b_1 \nu_h \frac{I_1^* E_h}{I_1} + \sum_{i=2}^{n-1} \left( b_i(\mu_h + \eta_i + \gamma_i)I_i^* - b_i \gamma_{i-1} \frac{I_i^* I_{i-1}}{I_i} \right) \\
& + b_n(\mu_h + \eta_n)I_n^* - b_n \gamma_{n-1} \frac{I_n^* I_{n-1}}{I_n}.
\end{aligned}
$$

Taking into account the construction of $b_i$, we have

$$ 2b_0(\mu_h + \nu_h)E_h^* = 2v_1(\mu_v + \delta_v)I_v^* = 2v_1 S_v^* \sum_{i=1}^n \bar{\beta}_i I_i^*, $$

$$ 
\begin{aligned}
b_i &= v_1 S_v^* \frac{\sum_{j=i}^n \bar{\beta}_j I_j^*}{(\mu_h + \eta_i + \gamma_i) I_i^*} \\
&= v_1 S_v^* \frac{\sum_{j=i}^n \bar{\beta}_j I_j^*}{\alpha_i I_i^*}, \quad \forall i = 1, 2, \dots, n-1,
\end{aligned}
$$

and,

$$ b_n = v_1 \bar{\beta}_n S_v^*. $$

It follows that:

$$ 
\begin{aligned}
\frac{\dot{W}}{v_1} ={}& \mathcal{A}_v + (b_0 \mathcal{A}_h / v_1) + 2S_v^* \sum_{i=1}^n \bar{\beta}_i I_i^* - S_v^* \sum_{i=1}^n \bar{\beta}_i I_i^* \frac{S_v^*}{S_v} - \frac{I_v^*}{I_v} S_v \sum_{i=1}^n \bar{\beta}_i I_i + 2S_v^* \sum_{i=1}^n \bar{\beta}_i I_i^* \\
& - S_v^* \sum_{i=1}^n \bar{\beta}_i I_i^* \frac{S_h^*}{S_h} - S_v^* \frac{\sum_{i=1}^n \bar{\beta}_i I_i^*}{S_h^* I_v^*} S_h \frac{E_h^* I_v}{E_h} + S_v^* \sum_{i=1}^n \bar{\beta}_i I_i^* - S_v^* \frac{\sum_{i=1}^n \bar{\beta}_i I_i^*}{E_h^*} \frac{I_1^* E_h}{I_1} \\
& + \sum_{i=2}^{n-1} \left( S_v^* \sum_{j=i}^n \bar{\beta}_j I_j^* - S_v^* \frac{\sum_{j=i}^n \bar{\beta}_j I_j^*}{I_{i-1}^*} \frac{I_i^* I_{i-1}}{I_i} \right) + \bar{\beta}_n S_v^* I_n^* - \bar{\beta}_n S_v^* \frac{I_n^*}{I_{n-1}^*} \frac{I_n^* I_{n-1}}{I_n}.
\end{aligned}
$$
By separating the first index from the sum, we obtain:

$$
\begin{aligned}
\dot{W}/v_1 ={}& A_v + (b_0 A_h / v_1) + 5 S_v^* \bar{\beta}_1 I_1^* - S_v^* \bar{\beta}_1 I_1^* \frac{S_v^*}{S_v} - \frac{I_v^*}{I_v} S_v \bar{\beta}_1 I_1 \\
& - S_v^* \bar{\beta}_1 I_1^* \frac{S_h^*}{S_h} - S_v^* \bar{\beta}_1 I_1^* \frac{S_h E_h^* I_v}{S_h^* I_v^* E_h} - S_v^* \bar{\beta}_1 I_1^* \frac{I_1^* E_h}{E_h^* I_1} \\
& + 5 S_v^* \sum_{i=2}^n \bar{\beta}_i I_i^* - S_v^* \sum_{i=2}^n \bar{\beta}_i I_i^* \frac{S_v^*}{S_v} - \frac{I_v^*}{I_v} S_v \sum_{i=2}^n \bar{\beta}_i I_i \\
& - S_v^* \sum_{i=2}^n \bar{\beta}_i I_i^* \frac{S_h^*}{S_h} - S_v^* \sum_{i=2}^n \bar{\beta}_i I_i^* \frac{S_h E_h^* I_v}{S_h^* I_v^* E_h} - S_v^* \sum_{i=2}^n \bar{\beta}_i I_i^* \frac{I_1^* E_h}{E_h^* I_1} \\
& + \sum_{i=2}^n \left( S_v^* \sum_{j=i}^n \bar{\beta}_j I_j^* - S_v^* \sum_{j=i}^n \bar{\beta}_j I_j^* \frac{I_i^* I_{i-1}}{I_i I_{i-1}^*} \right).
\end{aligned}
$$

Therefore,

$$
\begin{aligned}
\dot{W}/v_1 ={}& A_v + (b_0 A_h / v_1) + S_v^* \bar{\beta}_1 I_1^* \left( 5 - \frac{S_v^*}{S_v} - \frac{S_v}{S_v^*} \frac{I_1}{I_1^*} \frac{I_v^*}{I_v} - \frac{S_h^*}{S_h} - \frac{S_h I_v E_h^*}{S_h^* I_v^* E_h} - \frac{E_h I_1^*}{E_h^* I_1} \right) \\
& + S_v^* \sum_{i=2}^n \bar{\beta}_i I_i^* \left( 5 - \frac{S_v^*}{S_v} - \frac{S_v I_i I_v^*}{S_v^* I_i^* I_v} - \frac{S_h^*}{S_h} - \frac{S_h I_v E_h^*}{S_h^* I_v^* E_h} - \frac{E_h I_1^*}{E_h^* I_1} \right) \\
& + \sum_{i=2}^n \left( S_v^* \sum_{j=i}^n \bar{\beta}_j I_j^* - S_v^* \sum_{j=i}^n \bar{\beta}_j I_j^* \frac{I_i^* I_{i-1}}{I_i I_{i-1}^*} \right).
\end{aligned}
$$

By using the properties of nested sums, $\sum_{i=2}^{n} \left( \sum_{j=i}^{n} u_j w_i \right) = \sum_{i=2}^{n} u_i \sum_{j=2}^{i} w_j$, we obtain,

$$
\sum_{i=2}^{n} \left( \sum_{j=i}^{n} \bar{\beta}_j I_j^* \right) = \sum_{i=2}^{n} \bar{\beta}_i I_i^* (i-1),
$$

and

$$
\sum_{i=2}^{n} \left( \sum_{j=i}^{n} \bar{\beta}_j I_j^* \frac{I_i^* I_{i-1}}{I_i I_{i-1}^*} \right) = \sum_{i=2}^{n} \bar{\beta}_i I_i^* \left( \sum_{i=2}^{i} \frac{I_i^* I_{i-1}}{I_i I_{i-1}^*} \right).
$$

Finally,

$$
\begin{aligned}
\dot{W}/v_1 ={}& A_v + (b_0 A_h / v_1) + S_v^* \bar{\beta}_1 I_1^* \left( 5 - \frac{S_v^*}{S_v} - \frac{S_v}{S_v^*} \frac{I_1}{I_1^*} \frac{I_v^*}{I_v} - \frac{S_h^*}{S_h} - \frac{S_h I_v E_h^*}{S_h^* I_v^* E_h} - \frac{E_h I_1^*}{E_h^* I_1} \right) \\
& + S_v^* \sum_{i=2}^n \bar{\beta}_i I_i^* \left( 4 + i - \frac{S_v^*}{S_v} - \frac{S_v I_i I_v^*}{S_v^* I_i^* I_v} - \frac{S_h^*}{S_h} - \frac{S_h I_v E_h^*}{S_h^* I_v^* E_h} - \frac{E_h I_1^*}{E_h^* I_1} - \sum_{i=2}^i \frac{I_i^* I_{i-1}}{I_i I_{i-1}^*} \right).
\end{aligned}
$$

By setting $v_1 = 1$ and replacing $A_v$ and $A_h$ by their respective values, the final expression of $\dot{W}$ is:

$$
\begin{aligned}
\dot{W} ={}& (\mu_v + \delta_v) S_v^* \left( 2 - \frac{S_v^*}{S_v} - \frac{S_v}{S_v^*} \right) + b_0 \mu_h S_h^* \left( 2 - \frac{S_h^*}{S_h} - \frac{S_h}{S_h^*} \right) \\
& + S_v^* \bar{\beta}_1 I_1^* \left( 5 - \frac{S_v^*}{S_v} - \frac{S_v}{S_v^*} \frac{I_1}{I_1^*} \frac{I_v^*}{I_v} - \frac{S_h^*}{S_h} - \frac{S_h I_v E_h^*}{S_h^* I_v^* E_h} - \frac{E_h I_1^*}{E_h^* I_1} \right) \\
& + S_v^* \sum_{i=2}^n \bar{\beta}_i I_i^* \left( 4 + i - \frac{S_v^*}{S_v} - \frac{S_v I_i I_v^*}{S_v^* I_i^* I_v} - \frac{S_h^*}{S_h} - \frac{S_h I_v E_h^*}{S_h^* I_v^* E_h} - \frac{E_h I_1^*}{E_h^* I_1} - \sum_{i=2}^i \frac{I_i^* I_{i-1}}{I_i I_{i-1}^*} \right)
\end{aligned}
\tag{1.17}
$$
with $b_0 = S_v^* \sum_{i=1}^n \bar{\beta}_i I_i^*$. By using the geometric-arithmetic mean inequality, we conclude that $\dot{W}$ is negative definite, which proves the global asymptotic stability of the endemic equilibrium. $\square$

This result of global stability of the endemic equilibrium is new and extend the results of [37] in which the endemic equilibrium is proven to be locally asymptotically stable if $R_0^2 > 1$ and $R_0^2$ 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_0^2$ acts as a sharp threshold for the disease in the sense that the latter dies out if $R_0^2 \le 1$ and persists whenever $R_0^2 > 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.

## 2 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 [33]. 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 [13, 14], the authors investigated the dynamics of Chagas' disease and WNV, respectively, with two host species. In modeling TBRF, Johnson et al. [27] 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

$$ \dot{S}_j = \Lambda_j - b_j(N_j, N_v)\beta_{vh}S_j \frac{I_v}{N_v} - \mu_h S_h, \quad (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:

$$ \dot{S}_v = \Lambda_v - \sum_{j=1}^{m} a_j \beta_{hv} S_v \frac{\sum_{i=1}^{n} \beta_{j,i} I_{j,i}}{N_j} - (\mu_v + \delta_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) = \frac{a_j N_v}{N_j}. $$
Hence, for $j = 1, 2, \dots, m$, Equation (2.1) leads to,

$$ \dot{S}_j = \Lambda_j - a_j \beta_{vh} S_j \frac{I_v}{N_j} - \mu_j S_j. $$

Therefore, the dynamics of a model that captures the interactions between *m* host species and an arthropod vector is given by:

$$
\left\{
\begin{array}{l}
\text{For } j = 1, 2, \dots, m : \\
\dot{S}_j = \Lambda_j - a_j \beta_{v,j} S_j \frac{I_v}{N_j} - \mu_j S_j \\
\dot{E}_j = a_j \beta_{v,j} S_j \frac{I_v}{N_j} - (\mu_j + \nu_j) E_{h,j} \\
\dot{I}_{j,1} = \nu_j E_{h,j} - (\mu_j + \eta_{j,1} + \gamma_{j,1}) I_{j,1} \\
\dot{I}_{j,2} = \gamma_{j,1} I_{j,1} - (\mu_j + \eta_{j,2} + \gamma_{j,2}) I_{j,2} \\
\vdots \\
\dot{I}_{j,n-1} = \gamma_{j,n-2} I_{j,n-2} - (\mu_j + \eta_{j,n-1} + \gamma_{j,n-1}) I_{j,n-1} \\
\dot{I}_{j,n} = \gamma_{j,n-1} I_{j,n-1} - (\mu_j + \eta_{j,n} + \gamma_{j,n}) I_{j,n} \\
\dot{S}_v = \Lambda_v - \sum_{j=1}^{m} a_j S_v \frac{\sum_{i=1}^{n} \beta_{j,i} I_{j,i}}{N_j} - (\mu_v + \delta_v) S_v \\
\dot{I}_v = \sum_{j=1}^{m} a_j S_v \frac{\sum_{i=1}^{n} \beta_{j,i} I_{j,i}}{N_j} - (\mu_v + \delta_v) I_v.
\end{array}
\right.
\qquad (2.2)
$$

Parameters in System (2.2) are described in Table 2.

<table><caption>Table 2: Description of the parameters used in System (2.2).</caption><thead><tr><th>Parameters</th><th>Description</th></tr></thead><tbody><tr><td>Λ<sub>j</sub></td><td>Recruitment of Host j</td></tr><tr><td>Λ<sub>v</sub></td><td>Recruitment of vectors</td></tr><tr><td>a<sub>j</sub></td><td>Bitting rate on Host j</td></tr><tr><td>μ<sub>j</sub></td><td>Host j's death rate</td></tr><tr><td>β<sub>v,j</sub></td><td>Infectiousness of Host j to mosquitoes per biting</td></tr><tr><td>β<sub>j,i</sub></td><td>Vector's Infectiousness to host j at stage i per biting</td></tr><tr><td>ν<sub>j</sub></td><td>Host j's rate at which the exposed individuals become infected</td></tr><tr><td>η<sub>j,i</sub></td><td>Per capita recovery rate for Host j at stage i</td></tr><tr><td>γ<sub>j,i</sub></td><td>Host j's per capita progression rate from stage i to i + 1</td></tr></tbody></table>

**Remark 2.1.** Some authors [7, 16, 44] 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) = \frac{a N_v}{\sum_{j=1}^{m} 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 ([27, 44]). 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 $\beta_{k,i} = 0$ for all $k \in K$ and $i = 1, 2, \dots, 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 $\beta_j = (\beta_{j,1}, \dots, \beta_{j,m})^T$, System (2.2) could be written in a more compact vectorial form as:

$$ 
\left\{
\begin{aligned}
\dot{S}_1 &= \Lambda_1 e_1 - a_1 \beta_{v,1} S_1 \frac{I_v}{\bar{N}_1} - \mu_1 S_1 \\
\dot{S}_2 &= \Lambda_2 e_1 - a_2 \beta_{v,2} S_2 \frac{I_v}{\bar{N}_2} - \mu_2 S_2 \\
&\vdots \\
\dot{S}_m &= \Lambda_m e_1 - a_m \beta_{v,m} S_m \frac{I_v}{\bar{N}_m} - \mu_m S_m \\
\dot{I}_1 &= a_1 \beta_{v,1} S_1 \frac{I_v}{\bar{N}_1} e_1 + A_1 \mathbf{I}_1 \\
\dot{I}_2 &= a_2 \beta_{v,2} S_2 \frac{I_v}{\bar{N}_2} e_1 + A_2 \mathbf{I}_2 \\
&\vdots \\
\dot{I}_m &= a_m \beta_{v,m} S_m \frac{I_v}{\bar{N}_m} e_1 + A_m \mathbf{I}_m \\
\dot{I}_v &= \sum_{j=1}^{m} a_j \frac{\bar{N}_v - I_v}{\bar{N}_{h,j}} \left\langle \begin{pmatrix} 0 \\ \beta_j \end{pmatrix} \middle| \mathbf{I}_j \right\rangle - (\mu_v + \delta_v) I_v,
\end{aligned}
\right.
\quad (2.3) $$

where $\mathbf{I}_j = (E_j, I_{j,1}, I_{j,2}, \dots, I_{j,n})^T$ is the vector of infected of Host j ($j = 1, 2, \dots, m$) and

$$ A_j = \begin{pmatrix} -(\nu_j + \mu_j) & 0 & \dots & \dots & 0 \\ \nu_j & -\alpha_{j,1} & \dots & \dots & 0 \\ 0 & \gamma_{j,1} & -\alpha_{j,2} & \dots & 0 \\ \vdots & \ddots & \ddots & \vdots & \vdots \\ 0 & \dots & \ddots & \gamma_{j,n-1} & -\alpha_{j,n} \end{pmatrix} $$

where $\alpha_{j,i} = \gamma_{j,i} + \eta_{j,i} + \mu_j$, for $j = 1, \dots, m$ and $i = 1, 2, \dots, n$. The term $\alpha_{j,i}$ represents the mean period for which infected population of Host j in stage i leave this stage.

## 2.1 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

$$ \Omega = \left\{ (S_j, E_j, I_{j,i}, I_v) \in \mathbb{R}_+^{m(n+2)+1} \mid S_j + E_j + \sum_{i=1}^{n} I_{j,i} \le \frac{\Lambda_j}{\mu_j}, \text{ for } j = 1, \dots, m, I_v \le \frac{\Lambda_v}{\mu_v} \right\} $$

is a compact attracting positively invariant set for system (2.3).
The disease-free equilibrium of Model (2.3) is $(S^*, 0_{\mathbb{R}^{(n+1)m+1}})$ with $S^* = [S_1^*, S_2^*, \dots, S_m^*] = [\frac{\Lambda_1}{\mu_1}, \frac{\Lambda_2}{\mu_2}, \dots, \frac{\Lambda_m}{\mu_m}]$. The vector field of Model (2.3) could be subdivided into new infections vector $\mathcal{F}$ and transition vector $\mathcal{V}$ where

$$ \mathcal{F}(\mathbf{I}_1, \dots, \mathbf{I}_m, I_v) = \begin{pmatrix} a_1 \beta_{v,1} S_1 \frac{I_v}{\bar{N}_1} e_1 \\ a_2 \beta_{v,2} S_2 \frac{I_v}{\bar{N}_2} e_1 \\ \vdots \\ a_m \beta_{v,m} S_m \frac{I_v}{\bar{N}_m} e_1 \\ \sum_{j=1}^m a_j \frac{\bar{N}_v - I_v}{\bar{N}_{h,j}} \langle \begin{pmatrix} 0 \\ \beta_j \end{pmatrix} | \mathbf{I}_j \rangle \end{pmatrix}, \mathcal{V}(\mathbf{I}_1, \dots, \mathbf{I}_m, I_v) = \begin{pmatrix} A_1 \mathbf{I}_1 \\ A_2 \mathbf{I}_2 \\ \vdots \\ A_m \mathbf{I}_m \\ -(\mu_v + \delta_v) I_v \end{pmatrix}. $$

The Jacobian matrices at the DFE of $\mathcal{F}$ and $\mathcal{V}$ are $F = D\mathcal{F}(\mathbf{I}_1, \mathbf{I}_2, \dots, \mathbf{I}_m, I_v)|_{\text{DFE}}$ and $V = D\mathcal{V}(\mathbf{I}_1, \mathbf{I}_2, \dots, \mathbf{I}_m, I_v)|_{\text{DFE}}$ and are given by:

$$ F = \begin{pmatrix} \mathbf{0}_{n+1,n+1} & \mathbf{0}_{n+1,n+1} & \cdots & \cdots & a_1 \beta_{v,1} e_1 \\ \mathbf{0}_{n+1,n+1} & \mathbf{0}_{n+1,n+1} & \cdots & \cdots & a_2 \beta_{v,2} e_1 \\ \vdots & \vdots & \vdots & \vdots & \vdots \\ \mathbf{0}_{n+1,n+1} & \mathbf{0}_{n+1,n+1} & \cdots & \cdots & a_m \beta_{v,m} e_1 \\ a_1 \frac{\bar{N}_v}{\bar{N}_1} (0, \boldsymbol{\beta}_1)^T & a_2 \frac{\bar{N}_v}{\bar{N}_2} (0, \boldsymbol{\beta}_2)^T & \cdots & a_m \frac{\bar{N}_v}{\bar{N}_m} (0, \boldsymbol{\beta}_m)^T & 0 \end{pmatrix}, $$

and

$$ V = \operatorname{diag}(A_1, A_2, \dots, A_m, -\mu_v - \delta_v). $$

It follows that

$$ -V^{-1} = \operatorname{diag} \left( -A_1^{-1}, -A_2^{-1}, \dots, -A_m^{-1}, \frac{1}{\mu_v + \delta_v} \right). $$

Thus, the next generation matrix is:

$$ -FV^{-1} = \begin{pmatrix} \mathbf{0}_{n+1,n+1} & \cdots & \cdots & \frac{a_1 \beta_{v,1}}{\mu_v + \delta_v} e_1 \\ \mathbf{0}_{n+1,n+1} & \cdots & \cdots & \frac{a_2 \beta_{v,2}}{\mu_v + \delta_v} e_1 \\ \vdots & \vdots & \vdots & \vdots \\ \mathbf{0}_{n+1,n+1} & \cdots & \cdots & \frac{a_m \beta_{v,m}}{\mu_v + \delta_v} e_1 \\ a_1 \frac{\bar{N}_v}{\bar{N}_1} (0, \boldsymbol{\beta}_1)^T (-A_1^{-1}) & \cdots & a_m \frac{\bar{N}_v}{\bar{N}_m} (0, \boldsymbol{\beta}_m)^T (-A_m^{-1}) & 0 \end{pmatrix}. $$

The matrix $(-FV^{-1})^2$ is of rank one, and its largest eigenvalue is:

$$ \begin{aligned} \mathcal{R}_0^2(m) &= \sum_{j=1}^{m} \frac{a_j^2 \beta_{v,j}}{\mu_v + \delta_v} \frac{\bar{N}_v}{\bar{N}_{h,j}} (0, \boldsymbol{\beta}_j)^T (-A_j^{-1}) e_1 \\ &:= \sum_{j=1}^{m} \mathcal{R}_0^2, \end{aligned} $$
where $\mathcal{R}_{0,j}^2$ is the basic reproduction number when Host $j$ is the only coupled host to the vector. A more explicit expression of $\mathcal{R}_0^2$ could be obtained by computing $A_j^{-1}$. Indeed,

$$ -A_j^{-1} = \begin{pmatrix} \frac{1}{\nu_j + \mu_j} & 0 & \dots & \dots & 0 \\ \frac{\nu_j}{\alpha_{j,1}(\nu_j + \mu_j)} & \frac{1}{\alpha_{j,1}} & \dots & \dots & 0 \\ \frac{\gamma_{j,1}\nu}{\alpha_{j,1}\alpha_{j,2}(\nu_j + \mu_j)} & \frac{\gamma_{j,1}}{\alpha_{j,1}\alpha_{j,2}} & \frac{1}{\alpha_{j,2}} & \dots & 0 \\ \vdots & \ddots & \ddots & \ddots & \vdots \\ \frac{\gamma_{j,1}\dots\gamma_{j,n-1}\nu_j}{\alpha_{j,1}\alpha_{j,2}\dots\alpha_{j,n}(\nu_j + \mu_j)} & \dots & \dots & \frac{\gamma_{j,n-1}}{\alpha_{j,n-1}\alpha_{j,n}} & \frac{1}{\alpha_{j,n}} \end{pmatrix}. \quad (2.4) $$

Hence, we deduce that,

$$ \begin{aligned} (0, \beta_j)^T (-A_j^{-1}) e_1 &= \beta_{j,1} \frac{\nu_j}{\alpha_{j,1}(\nu_j + \mu_j)} + \beta_{j,2} \frac{\gamma_{j,1}\nu_j}{\alpha_{j,1}\alpha_{j,2}(\nu_j + \mu_j)} + \dots \\ & \quad + \beta_{j,n} \frac{\gamma_{j,1} \dots \gamma_{j,n-1} \nu_j}{\alpha_{j,1}\alpha_{j,2} \dots \alpha_{j,n}(\nu_j + \mu_j)} \\ &= \frac{\nu_j}{\nu_j + \mu_j} \left( \frac{\beta_{j,1}}{\alpha_{j,1}} + \frac{\beta_{j,2}\gamma_{j,1}}{\alpha_{j,1}\alpha_{j,2}} + \dots + \frac{\beta_{j,n}\gamma_{j,1} \dots \gamma_{j,n-1}}{\alpha_{j,1}\alpha_{j,2} \dots \alpha_{j,n}} \right) \\ &= \frac{\nu_j}{\nu_j + \mu_j} \sum_{i=1}^{n} \frac{\gamma_{j,1}\gamma_{j,2} \dots \gamma_{j,i-1}}{\alpha_{j,1}\alpha_{j,2} \dots \alpha_{j,i}} \beta_{j,i} \quad \text{with} \quad \gamma_{j,0} = 1. \end{aligned} $$

Finally, the basic reproduction number of a multi-host-vector model with $n$ stages is given by

$$ \mathcal{R}_0^2(m) = \sum_{j=1}^{m} \mathcal{R}_{0,j}^2, $$

where

$$ \mathcal{R}_{0,j}^2 = \frac{a_j^2 \beta_{v,j}}{\mu_v + \delta_v} \frac{\bar{N}_v}{\bar{N}_{h,j}} \frac{\nu_j}{\nu_j + \mu_j} \sum_{i=1}^{n} \frac{\gamma_{j,1}\gamma_{j,2}\dots\gamma_{j,i-1}}{\alpha_{j,1}\alpha_{j,2}\dots\alpha_{j,i}} \beta_{j,i}. $$

This result generalizes the results of [13, 14, 27], for which $j = 2$, and [37] for which $j = 1$. Similar remarks hold for [7].

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 [16]. Indeed, in [16] the authors found that the basic reproduction number in multi-host models, when not coupled with an arthropod vector, increases with host diversity for density-dependent transmission and decreases for frequency-dependent transmission. Our model consists of frequency-dependent transmission, but $\mathcal{R}_0^2(m)$ increases with host diversity. This contrast from [16] 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 [7, 16, 44], the authors define $b(\sum_{j=1}^m N_j, N_v)$ as the number of bites on all hosts. The latter will lead to a factor of $\sum_{j=1}^m 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 $\mathcal{R}_0^2(m)$ coincides with that of [16].
## 2.2 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_0^2(m)$.

**Theorem 2.1.**

The DFE is globally asymptotically stable for the multi-host Model (2.3) whenever $R_0^2(m) \le 1$.

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

$$ \mathcal{V} = \sum_{j=1}^{m} p_j \mathcal{L}_j + I_v $$

where $\mathcal{L}_j = c_{j,0}E_j + c_{j,1}I_{j,1} + c_{j,2}I_{j,2} + \dots + c_{j,n}I_{j,n}$ and $p_j = a_j\beta_{j,v}\frac{\bar{N}_v}{\bar{N}_{h,j}}$. The coefficients $c_{j,i}$ for $i = 1, 2, \dots, n$ are similar to those chosen in the proof of Theorem 1.1. They are defined recursively as follows:

$$ \left\{ \begin{aligned} c_{j,0} &= \frac{\nu_j}{\nu_j + \mu_j} c_{j,1} \\ c_{j,1} &= \left( \frac{\beta_{j,1}}{\alpha_{j,1}} + \frac{\gamma_{j,1}}{\alpha_{j,1}\alpha_{j,2}} \beta_{j,2} + \dots + \frac{\gamma_{j,1}\dots\gamma_{j,n-1}}{\alpha_{j,1}\alpha_{j,2}\dots\alpha_{j,n}} \beta_{j,n} \right) \\ c_{j,i+1} &= \frac{1}{\gamma_{j,i}} (\alpha_{j,i} c_{j,i} - \beta_{j,i}), \quad \forall i = 1, \dots, n-2 \\ c_{j,n} &= \frac{\beta_{j,n}}{\alpha_{j,n}}. \end{aligned} \right. \qquad (2.5) $$

The derivative of $\mathcal{V}$ along the solutions of System (2.3) is:

$$ \dot{\mathcal{V}} = \sum_{j=1}^{m} p_j \dot{\mathcal{L}}_j + \dot{I}_v, $$

where:

$$ \dot{\mathcal{L}}_j = c_{j,0}\dot{E}_j + \sum_{i=1}^{n} c_{j,i}\dot{I}_{j,i}. $$

By using the induction relationship (2.5) and as in the proof of Theorem 1.1, we obtain

$$ \dot{\mathcal{L}}_j = c_{j,0}a_j\beta_{v,j}S_j\frac{I_v}{\bar{N}_j} - \sum_{i=1}^{n}\beta_{j,i}I_{j,i}. \qquad (2.6) $$

Thus, the derivative of the function $\mathcal{V}$ becomes:

$$ \begin{aligned} \dot{\mathcal{V}} &= \sum_{j=1}^{m} p_j \left( c_{j,0} a_j \beta_{v,j} S_j \frac{I_v}{\bar{N}_j} - \sum_{i=1}^{n} \beta_{j,i} I_{j,i} \right) + \left( \sum_{j=1}^{m} a_j (\bar{N}_v - I_v) \frac{\sum_{i=1}^{n} \beta_{j,i} I_{j,i}}{\bar{N}_j} - (\mu_v + \delta_v) I_v \right) \\ &= \sum_{j=1}^{m} \left( a_j \frac{\bar{N}_v}{\bar{N}_j} - p_j \right) \sum_{i=1}^{n} \beta_{j,i} I_{j,i} - \sum_{j=1}^{m} a_j I_v \frac{\sum_{i=1}^{n} \beta_{j,i} I_{j,i}}{\bar{N}_j} \\ &\quad + \left( \sum_{j=1}^{m} p_j c_{j,0} a_j \beta_{v,j} \frac{S_j}{\bar{N}_j} - (\mu_v + \delta_v) \right) I_v \\ &= - \sum_{j=1}^{m} a_j I_v \frac{\sum_{i=1}^{n} \beta_{j,i} I_{j,i}}{\bar{N}_j} + \left( \sum_{j=1}^{m} c_{j,0} a_j^2 \beta_{v,j} \frac{\bar{N}_v S_j}{\bar{N}_j \bar{N}_j} - (\mu_v + \delta_v) \right) I_v, \quad \text{since } p_j = a_j \frac{\bar{N}_v}{\bar{N}_j}. \end{aligned} $$
Hence, we obtain,

$$
\begin{align}
\dot{V} &= -\sum_{j=1}^{m} a_j \beta_{j,v} I_v \frac{\sum_{i=1}^{n} I_{j,i}}{\bar{N}_j} + (\mu_v + \delta_v) \left( \sum_{j=1}^{m} \mathcal{R}_{0,j}^2 \frac{S_j}{\bar{N}_j} - 1 \right) I_v \nonumber \\
&\leq -I_v \sum_{j=1}^{m} a_j \beta_{j,v} \frac{\sum_{i=1}^{n} I_{j,i}}{\bar{N}_j} + (\mu_v + \delta_v) (\mathcal{R}_0^2(m) - 1) I_v, \quad \text{since } S_j \leq \bar{N}_j, \forall j = 1, \dots, m, \tag{2.7} \\
&\leq 0 \quad \text{whenever } \mathcal{R}_0^2(m) \leq 1. \nonumber
\end{align}
$$

It is not difficult to see that the largest invariant set within {$\dot{V} = 0$} is reduced to the DFE, which is therefore globally asymptotically stable in the compact set $\Omega$ thanks to ([4], Theorem 3.7.11, page 346, [30] Theorem 3). Since $\Omega$ is an attractive set, it follows that the DFE is globally asymptotically stable in $\mathbb{R}^{m(n+2)+1}$. $\square$

**Theorem 2.2.**

A unique endemic equilibrium exists for the multi-host Model (2.3) whenever $\mathcal{R}_0^2(m) > 1$.

Proof.

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

$$E_j = \frac{\alpha_{j,1}}{\nu_j} I_{j,1}^*, \quad I_{j,1}^* = \frac{\alpha_{j,2}}{\gamma_{j,1}} I_{j,2}^*,$$

and for $i = 2, \dots, n$,

$$I_{j,i-1}^* = \frac{\alpha_{j,i}}{\gamma_{j,i-1}} I_{j,i}^*.$$

Hence, as in (1.11) for the one host case, we can show that:

$$\sum_{i=1}^{n} \beta_{j,i} I_{j,i}^{*} = \alpha_{j,1} c_{j,1} I_{j,1}^{*}. \qquad (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^* = \frac{\mu_j \bar{N}_j^2}{\mu_j \bar{N}_j + a_j \beta_{v,j} I_v^*}$$

$$
\begin{align*}
E_j^* &= \frac{a_j \beta_{v,j} \mu_j \bar{N}_j^2}{\mu_j + \nu_j \mu_j \bar{N}_j + a_j \beta_{v,j} I_v^* \bar{N}_j} \frac{I_v^*}{\bar{N}_j} \\
&= \frac{a_j \beta_{v,j} \mu_j \bar{N}_j I_v^*}{\mu_j + \nu_j \mu_j \bar{N}_j + a_j \beta_{v,j} I_v^*}
\end{align*}
$$

$$
\begin{align*}
I_{j,1}^* &= \frac{\nu_j}{\alpha_{j,1}} E_j^* \\
&= \frac{\nu_j a_j \beta_{v,j} \mu_j \bar{N}_j I_v^*}{\alpha_{j,1} \mu_j + \nu_j \mu_j \bar{N}_j + a_j \beta_{v,j} I_v^*}.
\end{align*}
$$

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
$$
\begin{align*}
0 &= (\bar{N}_v - I_v^*) \sum_{j=1}^{m} a_j \frac{\sum_{i=1}^{n} \beta_{j,i} I_{j,i}^*}{\bar{N}_j} - (\mu_v + \delta_v) I_v^* \\
&= (\bar{N}_v - I_v^*) \sum_{j=1}^{m} a_j \frac{\alpha_{j,1} c_{j,1} I_{j,1}^*}{\bar{N}_j} - (\mu_v + \delta_v) I_v^* \\
&= (\bar{N}_v - I_v^*) \sum_{j=1}^{m} a_j \frac{\alpha_{j,1} c_{j,1} \nu_j}{\bar{N}_j \alpha_{j,1} \mu_j + \nu_j \mu_j \bar{N}_j + a_j \beta_{v,j} I_v^*} - (\mu_v + \delta_v) I_v^* \\
&= (\bar{N}_v - I_v^*) \sum_{j=1}^{m} a_j c_{j,1} \frac{\nu_j a_j \beta_{v,j}}{\mu_j + \nu_j \mu_j \bar{N}_j + a_j \beta_{v,j} I_v^*} - (\mu_v + \delta_v) I_v^* \tag{2.9}
\end{align*}
$$

Since the equilibrium is endemic $I_v^* > 0$, Equation (2.9) is satisfied if, and only if:

$$
\begin{align}
0 &= \frac{\bar{N}_v - I_v^*}{\mu_v + \delta_v} \sum_{j=1}^{m} c_{j,1} \frac{\nu_j a_j^2 \beta_{v,j}}{\mu_j + \nu_j \mu_j \bar{N}_j + a_j \beta_{v,j} I_v^*} - 1 \nonumber \\
  &:= f(I_v^*). \tag{2.10}
\end{align}
$$

The function $f(I_v^*)$ is decreasing and $f(0) = \mathcal{R}_0^2(m) - 1$. Moreover, $f(\bar{N}_v) = -1 < 0$. Hence, the equation $f(I_v^*)$ has a unique positive root in $(0, \bar{N}_v)$ if, and only if $f(0) > 0$, that is, whenever $\mathcal{R}_0^2(m) > 1$.

Since the other components $S_j^*, E_j^*, I_{i,j}^*$, for $j = 1, 2, \dots, m$ and $i = 1, 2, \dots, n$, are uniquely expressed in terms of $I_v^*$, we conclude that there is a unique endemic equilibrium in Int($\Omega$) whenever $\mathcal{R}_0^2(m) > 1$. $\square$

**Theorem 2.3.** The endemic equilibrium is globally asymptotically stable on $\mathbb{R}^{m(n+2)+1}$.

*Proof.*

We consider the following Lyapunov function

$$
\mathcal{V} = \sum_{j=1}^{m} A_{j0} (S_j - S_j^* \log S_j + E_j - E_j^* \log E_j) + \sum_{j=1}^{m} \left( \sum_{i=1}^{n} A_{j,i} (I_{j,i} - I_{j,i}^* \log I_{j,i}) \right) \\
+ V_1 (S_v - S_v^* \log S_v) + V_1 (I_v - I_v^* \log I_v) \quad (2.11)
$$

where the coefficients $V_1, A_{j0}, A_{ji}$, for $i = 1, 2, \dots, n$ and $j = 1, 2, \dots, m$ are positive to be determined later. The function $\mathcal{V}$ is thus positive-definite.

The derivative of $\mathcal{V}$ along trajectories is:

$$
\begin{align}
\dot{\mathcal{V}} ={}& -V_1(\mu_v + \delta_v)I_v + \mathcal{A}_v \nonumber \\
& + V_1 \sum_{j=1}^{m} \left( a_j \frac{S_v^*}{\bar{N}_j} \sum_{i=1}^{n} \beta_{j,i} I_{j,i}^* \left( 2 - \frac{S_v^*}{S_v} - \frac{S_v}{S_v^*} \frac{I_{j,i}}{I_{j,i}^*} \frac{I_v^*}{I_v} \right) + a_j S_v^* \frac{\sum_{i=1}^{n} \beta_{j,i} I_{j,i}}{\bar{N}_j} \right) \nonumber \\
& + \sum_{j=1}^{m} \left[ A_{hj} + A_{j0} a_j \beta_{v,j} S_j^* \frac{I_v^*}{\bar{N}_j} \left( 2 - \frac{S_j^*}{S_j} - \frac{S_j}{S_j^*} \frac{I_v}{I_v^*} \frac{E_j^*}{E_j} \right) + A_{j0} \left( a_j \beta_{v,j} S_j^* \frac{I_v}{\bar{N}_j} \right) \right. \nonumber \\
& \left. - A_{j0}(\mu_v + \nu_j)E_j + A_{j1} (\nu_j E_j - \alpha_{j,1} I_{j,1}) - A_{j1} \nu_j E_j \left( \frac{I_{j,1}^*}{I_{j,1}} \right) + A_{j1} \alpha_{j,1} I_{j,1}^* \right. \nonumber \\
& \left. + \sum_{i=2}^{n} \left( A_{ji} (\gamma_{j,i-1} I_{j,i-1} - \alpha_{j,i} I_{j,i}) - A_{ji} \gamma_{j,i-1} I_{j,i-1} \left( \frac{I_{j,i}^*}{I_{j,i}} \right) + A_{ji} \alpha_{j,i} I_{j,i}^* \right) \right], \tag{2.12}
\end{align}
$$
where $\mathcal{A}_v = (\mu_v + \delta_v)S_v^* \left(2 - \frac{S_v^*}{S_v} - \frac{S_v}{S_v^*}\right)$ and $\mathcal{A}_{hj} = \mu_h S_j^* \left(2 - \frac{S_j^*}{S_j} - \frac{S_j}{S_j^*}\right)$.

We choose the coefficients $A_{ji}$ in a similar fashion as in the one host case, that is,

$$ 
\begin{aligned}
A_{j0}(\mu_j + \nu_j)E_j^* &= A_{j0}a_j\beta_{v,j}S_j^*\frac{I_v^*}{N_j} \\
&= V_1a_j\frac{S_v^*}{N_j}\left(\sum_{i=1}^n\beta_{j,i}I_{j,i}^*\right) \quad \text{for all } j = 1, 2, \dots, m,
\end{aligned} 
\quad (2.13) $$

and

$$ A_{j1} = \frac{1}{\nu_j E_j^*} V_1 a_j \frac{S_v^*}{N_j} \left( \sum_{i=1}^n \beta_{j,i} I_{j,i}^* \right) \quad \text{for all } j = 1, 2, \dots, m. $$

Also,

$$ A_{ji} = \frac{\sum_{k=i}^{n} \beta_{j,k} I_{j,k}^*}{\alpha_{j,i} I_{j,i}^*} V_1 a_j \beta_{j,v} \frac{S_v^*}{N_j} \quad \text{for all } 2 \le i \le n-1, 1 \le j \le m, $$

and

$$ \alpha_{j,n} A_{jn} = V_1 a_j \beta_{j,v} \frac{S_v^*}{N_j} \quad \text{for all } 1 \le j \le m. $$

We obtain

$$ 
\begin{aligned}
\dot{\mathcal{V}} ={}& V_1 \sum_{j=1}^{m} a_j \beta_{j,1} \frac{S_v^*}{N_j} I_{j,1}^* \left( 5 - \frac{S_v^*}{S_v} - \frac{S_v}{S_v^*} \frac{I_{j,1}}{I_v} \frac{I_v^*}{I_v} - \frac{S_j^*}{S_j} - \frac{S_j}{S_j^*} \frac{I_v}{I_v^*} \frac{E_j^*}{E_j} - \frac{E_j}{E_j^*} \frac{I_{j,1}^*}{I_{j,1}} \right) \\
& + V_1 \sum_{j=1}^{m} a_j \frac{S_v^*}{N_j} \sum_{i=2}^{n} \beta_{j,i} I_{j,i}^* \left( 4 + i - \frac{S_v^*}{S_v} - \frac{S_v}{S_v^*} \frac{I_{j,i}}{I_{j,i}^*} \frac{I_v^*}{I_v} - \frac{S_j^*}{S_j} - \frac{S_j}{S_j^*} \frac{I_v}{I_v^*} \frac{E_j^*}{E_j} - \frac{E_j}{E_j^*} \frac{I_{j,1}^*}{I_{j,1}} \right. \\
& \qquad \left. - \sum_{k=2}^{i} \frac{I_{j,k-1}}{I_{j,k-1}^*} \frac{I_{j,k}^*}{I_{j,k}} \right) \\
& + \mathcal{A}_v + \sum_{j=1}^{m} \mathcal{A}_{hj}.
\end{aligned}
 $$

The function $\dot{\mathcal{V}}$ is therefore negative-definite. Since $\mathcal{V}$ is positive-definite, the global stability of the interior endemic equilibrium on $\mathbb{R}^{m(n+2)+1}$ follows from the Lyapunov stability theorem. $\square$

The result of global stability of the endemic equilibrium generalizes [13, 14] for an arbitrary number of hosts. Indeed, in [13], 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 an host's infectious stages. Similarly, [14] 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 $\mathcal{R}_0^2(m) \le 1$ and persists otherwise. However, $\mathcal{R}_0^2(m) = \sum_{j=1}^m \mathcal{R}_{0,j}^2$, where $\mathcal{R}_{0,j}$ is the basic reproduction number in presence of Host j only. Hence,

$$ \mathcal{R}_0^2(m) \le 1 \implies \mathcal{R}_{0,j}^2 \le 1 \text{ for all } j = 1, 2, \dots, m. $$
Thus, the eradication of 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.

## 2.3 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 $\mathcal{R}_0^2(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 considered separately, the system behave as two isolated host-vector systems, with the same parameter values, the corresponding basic reproduction numbers are $\mathcal{R}_{0,1}^2 = 0.66$ and $\mathcal{R}_{0,2}^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 and 3 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 $\beta$s lower than those given in Table 3 will lead to the overall basic reproduction number $\mathcal{R}_0^2(2) \le 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 $\mathcal{R}_{0,1}^2 > 1$ and $\mathcal{R}_{0,2}^2 > 1$, the trajectories of the isolated host-vector systems converge to an endemic equilibria as expected (Theorem 1.2). We decided not to display the figures of these cases.

# 3 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 $\mathcal{R}_0$. Our results show that the dynamics of the model is robust and is tied to $\mathcal{R}_0$. The disease-free equilibrium is globally asymptotically stable (GAS) if $\mathcal{R}_0 \le 1$ and unstable otherwise (Theorem 1.1). A unique endemic equilibrium exists and is GAS whenever $\mathcal{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 $\mathcal{R}_0^2(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, $\mathcal{R}_{0,j}^2$, for $j = 1, 2, ..., m$. We proved that the system with $m$ hosts remains robust as its asymptotic behavior is solely determined by $\mathcal{R}_0(m)$. The disease dies out
<figure><img src="image_3.png" /><figcaption>Figure 2: Simulation of two separate systems with one host each and four stages of infection. The $R_{0,j}^2$ 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.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>Figure 3: Simulation of a two-host system where there are four stages of infection for each host. The $R_{0,j}^2$ values are 0.66 for the host 1 population and 0.84 for the host 2 population, giving a combined system $R_0^2(2) = 1.50$. The infected populations eventually go to the nonzero endemic equilibrium.</figcaption></figure>

from all hosts if $R_0(m) \le 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)$
<table><caption>Table 3: Parameters used in the simulations.</caption><thead><tr><th>Parameter</th><th>Host 1</th><th>Host 2</th></tr></thead><tbody><tr><td>Λ<sub>j</sub></td><td>100</td><td>100</td></tr><tr><td>Λ<sub>v</sub></td><td>10000</td><td>10000</td></tr><tr><td>β<sub>v,j</sub></td><td>0.2</td><td>0.2</td></tr><tr><td>μ<sub>j</sub></td><td>1/28875</td><td>1/28875</td></tr><tr><td>μ<sub>j</sub></td><td>1/15</td><td>1/15</td></tr><tr><td>ν<sub>j</sub></td><td>1/15</td><td>1/15</td></tr><tr><td>a<sub>j</sub></td><td>0.37</td><td>0.37</td></tr><tr><td>δ<sub>v</sub></td><td>0.001</td><td>0.001</td></tr><tr><td>β<sub>j,1</sub></td><td>0.2</td><td>0.2</td></tr><tr><td>β<sub>j,2</sub></td><td>0.3</td><td>0.3</td></tr><tr><td>β<sub>j,3</sub></td><td>0.4</td><td>0.4</td></tr><tr><td>β<sub>j,4</sub></td><td>0.5</td><td>0.5</td></tr><tr><td>η<sub>j,1</sub></td><td>0.01</td><td>0.01</td></tr><tr><td>η<sub>j,2</sub></td><td>0.007</td><td>0.007</td></tr><tr><td>η<sub>j,3</sub></td><td>0.004</td><td>0.004</td></tr><tr><td>η<sub>j,4</sub></td><td>0.001</td><td>0.001</td></tr><tr><td>γ<sub>j,1</sub></td><td>0.125</td><td>0.125</td></tr><tr><td>γ<sub>j,2</sub></td><td>0.01</td><td>0.0875</td></tr><tr><td>γ<sub>j,3</sub></td><td>0.05</td><td>0.05</td></tr><tr><td>γ<sub>j,4</sub></td><td>0.025</td><td>0.0125</td></tr><tr><td>R<sup>2</sup><sub>0,j</sub></td><td>0.66</td><td>0.84</td></tr><tr><td>R<sup>2</sup><sub>0</sub>(2)</td><td></td><td>1.50</td></tr></tbody></table>

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 than unity even if $R_{0,j} < 1$ for $j = 1, 2, \dots, m$, that is, for all hosts. Therefore, host heterogeneity favors the spread of the infection, and thus causes the so-called “amplification effect” [34] and opposed to the “dilution effect” [34, 36, 43], for which the increase in host diversity could potentially decrease or drive out the pathogens [36]. As $R_0^2(m) = \sum_{j=1}^{m} R_{0,j}^2$ 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 [12] concept, for which human, animal and ecosystem health are ecologically interconnected [35].
References

[1] R. M. ANDERSON AND R. M. MAY, *Infectious Diseases of Humans*. *Dynamics and Control*, Oxford science publications, 1991.

[2] M. BEGON AND R. G. BOWERS, *Host-host-pathogen models and microbial pest control: the effect of host self regulation*, Journal of Theoretical Biology, 169 (1994), pp. 275-287.

[3] M. BEGON, R. G. BOWERS, N. KADIANAKIS, AND D. E. HODGKINSON, *Disease and community structure: the importance of host self-regulation in a host-host-pathogen model*, The American Naturalist, 139 (1992), pp. 1131-1150.

[4] N. P. BHATIA AND G. P. SZEGÖ, *Dynamical systems: Stability theory and applications*, Lecture Notes in Mathematics, No. 35, Springer Verlag, Berlin, 1967.

[5] N. P. BHATIA AND G. P. SZEGÖ, *Stability Theory of Dynamical Systems*, Springer-Verlag, 1970.

[6] R. G. BOWERS AND M. BEGON, *A host-host-pathogen model with free-living infective stages*, applicable to microbial pest control, Journal of Theoretical Biology, 148 (1991), pp. 305-329.

[7] C. BOWMAN, A. GUMEL, P. VAN DEN DRIESSCHE, J. WU, AND H. ZHU, *A mathematical model for assessing control strategies against West Nile virus*, Bull Math Biol, 67 (2005), pp. 1107-1133.

[8] C. CASTILLO-CHAVEZ AND H. R. THIEME, *Asymptotically autonomous epidemic models*, in Mathematical Population Dynamics: Analysis of Heterogeneity, Volume One: Theory of Epidemics, O. Arino, A. D.E., and M. Kimmel, eds., Wuerz, 1995.

[9] CENTERS FOR DISEASE CONTROL AND PREVENTION, *Tickborne diseases of the united states*, Centers for Disease Control and Prevention, https://www.cdc.gov/lyme/resources/tickbornediseases.pdf, (2015).

[10] ——, *Transmission of japanese encephalitis virus*, Centers for Disease Control and Prevention, https://www.cdc.gov/japaneseencephalitis/symptoms/index.html, (2015).

[11] ——, *Zoonotic diseases*, Centers for Disease Control and Prevention, https://www.cdc.gov/onehealth/basics/zoonotic-diseases.html, (2015).

[12] ——, *One health*, Centers for Disease Control and Prevention, https://www.cdc.gov/onehealth/, (2017).

[13] G. CRUZ-PACHECO, L. ESTEVA, AND C. VARGAS, *Control measures for chagas disease*, Mathematical biosciences, 237 (2012), pp. 49-60.

[14] ——, *Multi-species interactions in west nile virus infection*, Journal of biological dynamics, 6 (2012), pp. 281-298.

[15] O. DIEKMANN AND J. A. P. HEESTERBEEK, *Mathematical epidemiology of infectious diseases*, Wiley Series in Mathematical and Computational Biology, John Wiley & Sons Ltd., Chichester, 2000. Model building, analysis and interpretation.

[16] A. DOBSON, *Population dynamics of pathogens with multiple host species*, the american naturalist, 164 (2004), pp. S64-S78.

[17] M. S. DWORKIN, T. G. SCHWAN, D. E. ANDERSON, AND S. M. BORCHARDT, *Tick-borne relapsing fever*, Infectious disease clinics of North America, 22 (2008), pp. 449-468.

[18] J. R. FRANCO, P. P. SIMARRO, A. DIARRA, AND J. G. JANNIN, *Epidemiology of human african trypanosomiasis*, Clin Epidemiol, 6 (2014), pp. 257-275.

[19] J. GREENMAN AND P. HUDSON, *Parasite-mediated and direct competition in a two-host shared macroparasite system*, Theoretical Population Biology, 57 (2000), pp. 13-34.
[20] H. GUO AND M. LI, *Global dynamics of a staged progression model for infectious diseases*, Math. Biosci. Eng., 3 (2006), pp. 513-525.

[21] H. GUO AND M. Y. LI, *Global dynamics of a staged progression model for infectious diseases*, Mathematical Biosciences and Engineering, 3 (2006), p. 513.

[22] H. GUO, M. Y. LI, AND Z. SHUAI, *Global dynamics of a general class of multistage models for infectious diseases*, SIAM Journal on Applied Mathematics, 72 (2012), pp. 261-279.

[23] R. D. HOLT AND J. PICKERING, *Infectious disease and species coexistence: a model of lotka-volterra form*, The American Naturalist, 126 (1985), pp. 196-211.

[24] M. J. HOMER, I. AGUILAR-DELFIN, S. R. TELFORD, P. J. KRAUSE, AND D. H. PERSING, *Babesiosis*, Clinical microbiology reviews, 13 (2000), pp. 451-469.

[25] A. IGGIDR, J. MBANG, G. SALLET, AND J. TEWA, *multi-compartment models*, DCDS, supplement (2007), pp. 506-519.

[26] C. K. JOHNSON, P. L. HITCHENS, T. S. EVANS, T. GOLDSTEIN, K. THOMAS, A. CLEMENTS, D. O. JOLY, N. D. WOLFE, P. DASZAK, W. B. KARESH, AND J. K. MAZET, *Spillover and pandemic properties of zoonotic viruses with high host plasticity*, Scientific reports, 5 (2015), p. 14830.

[27] T. L. JOHNSON, E. L. LANDGUTH, AND E. F. STONE, *Modeling relapsing disease dynamics in a host-vector community*, PLoS Negl Trop Dis, 10 (2016), p. e0004428.

[28] W. O. KERMACK AND A. G. MCKENDRICK, *Contributions to the mathematical theory of epidemics-i.*, 1927., Bull Math Biol, 53 (1991), pp. 33-55.

[29] S. U. KHAN, H. SALJE, A. HANNAN, M. A. ISLAM, A. M. BHUYAN, M. A. ISLAM, M. Z. RAHMAN, N. NAHAR, M. J. HOSSAIN, S. P. LUBY, ET AL., *Dynamics of japanese encephalitis virus transmission among pigs in northwest bangladesh and the potential impact of pig vaccination*, PLoS Negl Trop Dis, 8 (2014), p. e3166.

[30] J. LASALLE, *Stability theory for ordinary differential equations*, J. Differ. Equations, 4 (1968), pp. 57-65.

[31] J. P. LASALLE AND S. LEFSCHETZ, *Stability by Liapunov's direct method*, Academic Press, 1961.

[32] J. O. LLOYD-SMITH, D. GEORGE, K. M. PEPIN, V. E. PITZER, J. R. PULLIAM, A. P. DOBSON, P. J. HUDSON, AND B. T. GRENFELL, *Epidemic dynamics at the human-animal interface*, science, 326 (2009), pp. 1362-1367.

[33] A. A. MARFIN, L. R. PETERSEN, M. EIDSON, J. MILLER, J. HADLER, C. FARELLO, B. WERNER, G. L. CAMPBELL, M. LAYTON, P. SMITH, ET AL., *Widespread west nile virus activity, eastern united states, 2000.*, Emerging infectious diseases, 7 (2001), p. 730.

[34] E. MILLER AND A. HUPPERT, *The effects of host diversity on vector-borne disease: the conditions under which diversity will amplify or dilute the disease risk*, PLoS One, 8 (2013), p. e80279.

[35] ONE HEALTH, *One health*, One Health Commission, https://www.onehealthcommission.org/, (2017).

[36] R. S. OSTFELD AND F. KEESING, *Effects of host diversity on infectious disease*, Annual Review of Ecology, Evolution, and Systematics, 43 (2012), pp. 157-182.

[37] C. PALMER, E. LANDGUTH, E. STONE, AND T. JOHNSON, *The dynamics of vector-borne relapsing diseases*, arXiv preprint arXiv:1605.07647, (2016).

[38] T. R. PARK, *Age-dependence in Epidemic Models of Vector-borne Infections*, PhD thesis, University of Alabama in Huntsville, 2004.

[39] R. ROSS, *The prevention of malaria*, John Murray, 1911.
[40] R. ROSS, *An application of the theory of probabilities to the study of a priori pathometry*. part i, Proc. R. Soc. Lond. A, 92 (1916), pp. 204-230.

[41] R. ROSS AND H. HUDSON, *An application of the theory of probabilities to the study of a priori pathometry*. part ii, Proc. R. Soc. Lond. A, 93 (1917), pp. 212-225.

[42] R. ROSS AND H. HUDSON, *An application of the theory of probabilities to the study of a priori pathometry*. part iii, Proc. R. Soc. Lond. A, 93 (1917), pp. 225-240.

[43] D. J. SALKELD, K. A. PADGETT, AND J. H. JONES, *A meta-analysis suggesting that the relationship between biodiversity and risk of zoonotic pathogen transmission is idiosyncratic*, Ecology Letters, 16 (2013), pp. 679-686.

[44] J. E. SIMPSON, P. J. HURTADO, J. MEDLOCK, G. MOLAEI, T. G. ANDREADIS, A. P. GALVANI, AND M. A. DIUK-WASSER, *Vector host-feeding preferences drive transmission of multi-host pathogens: West nile virus as a model system*, Proceedings of the Royal Society of London B: Biological Sciences, (2011), p. rspb20111282.

[45] P. M. SOUTHERN JR AND J. P. SANFORD, *Relapsing fever: a clinical and microbiological review*, Medicine, 48 (1969), pp. 129-150.

[46] L. H. TAYLOR, S. M. LATHAM, AND E. MARK, *Risk factors for human disease emergence*, Philosophical Transactions of the Royal Society of London B: Biological Sciences, 356 (2001), pp. 983-989.

[47] P. VAN DEN DRIESSCHE AND J. WATMOUGH, *reproduction numbers and sub-threshold endemic equilibria for compartmental models of disease transmission*, Math. Biosci., (2002), pp. 29-48.

[48] J. X. VELASCOHERNANDEZ, *A model for chagas disease involving transmission by vectors and blood transfusion*, Theoretical population biology, 46 (1994), pp. 1-31.

[49] M. VIDYASAGAR, *Decomposition techniques for large-scale systems with nonadditive interactions: Stability and stabilizability*, IEEE Trans. Autom. Control, 25 (1980), pp. 773-779.

[50] O. WENGERT, M. KOPP, E. SIEBERT, W. STENZEL, G. HEGASY, N. SUTTORP, A. STICH, AND T. ZOLLER, *Human african trypanosomiasis with 7-year incubation period: clinical, laboratory and neuroimaging findings*, Parasitology international, 63 (2014), pp. 557-560.

[51] WORLD BANK, *People, pathogens and our planet, volume 1: Towards a one health approach for controlling zoonotic diseases*, World Bank, Agriculture and Rural Development Health, Nutrition and Population, (2010).