Global analysis of multi-host and multi-vector epidemic models.

We formulate a multi-group and multi-vector epidemic model in which hosts' dynamics is captured by staged-progression S E I R framework and the dynamics of vectors is captured by an SI framework. The proposed model describes the evolution of a class of zoonotic infections where the pathogen is shared by m host species and transmitted by p arthropod vector species. In each host, the infectious period is structured into n stages with a corresponding infectiousness parameter to each vector species. We determine the basic reproduction number R 0 2 ( m , n , p ) and investigate the dynamics of the systems when this threshold is less or greater than one. We show that the dynamics of the multi-host, multi-stage, and multi-vector system is completely determined by the basic reproduction number and the structure of the host-vector network configuration. Particularly, we prove that the disease-free equilibrium is globally asymptotically stable (GAS) whenever R 0 2 ( m , n , p ) < 1 , and a unique strongly endemic equilibrium exists and is GAS if R 0 2 ( m , n , p ) > 1 and the host-vector configuration is irreducible. That is, either the disease dies out or persists in all hosts and all vector species.

Introduction

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 [27], [41], 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 [45]. Furthermore, it is estimated that zoonoses cause over 20 billion and 200 billion USD of direct and indirect economic burden across the world respectively [27], [45]. Although the morbidity and mortality of most ID decreased, the incidence of zoonoses have increased [28]. 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 [25]. 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], [30]. Recent research has been focused on understanding the dynamics and control of vector-borne zoonoses with one vector and two hosts (see [35], [39] for WNV, [14] for Chagas' disease, and [26], [33] 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 transmitting 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 [20], [34]. Another example of a highly complex zoonosis is the eastern equine encephalitis virus (EEEV). Indeed, as illustrated in [13], [31], 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 [31], 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 [29]. 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 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 . 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 [38], [24] and sexually transmitted diseases [22]. The global dynamics of vector-borne diseases with one host species and one vector species have been investigated in different settings (see [10] for a model with time delays, [11] for a competitive exclusion result on a multi-strain model with distributed delay and [9] for the interplay between sterile and wild mosquitoes populations). For our system, a multi-host, multi-stage and multi-vector model, we will prove indeed that it has a unique strongly EE (in the sense of Thieme [42]) under the assumption that the host-vector network configuration is irreducible and . 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 [17], [18] 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 [36], [37], the authors used these tools to investigate for stage-progression models for water-borne diseases and the authors in [23] successfully generalized the approach to a multi-group vector-borne 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.

Formulation of the model

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 structure where n designates the number of infectious stages in the evolution of the disease within each host species. Let represent the total population of each host species i () and let , , , and denote the susceptible, latent, infectious at stage l (), and recovered populations of host species i, respectively. Let be the total population of arthropod vectors of species j, (), each of which is composed by susceptible vectors and infectious vectors . The susceptible populations of host species i are generated through a constant recruitment , subjected to a natural death rate of 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 . 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 , progressing to the next stage of infection at a rate or naturally die at a rate . The susceptible arthropod vectors of species j, , could be infected by all infectious hosts, of any species and of any stage, at different rates. Moreover, they are replenished at a constant rate and die either by natural death, at rate , or due to control strategies specific for each vector of species j, at a rate .

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 Fig. 1 and the parameters are described in Table 1 . The total population of host species and vector species are asymptotically constant as their dynamics are given by and , respectively. Hence, by using the theory of asymptotic systems [12], [44] and by denoting the limits of and by and , System (1) is asymptotically equivalent to its limit system. Moreover, in this case, it could be written in a compact form as: where , , , is the vector of infectious at stage k () for all hosts species. Also, and for all i. For vectors, and denote the vectors of susceptible and infected, respectively. The matrices A, and are given by for . A represents the biting/landing rates matrix, captures the infectiousness of host species to vector species; and the matrices represents the infectiousness of vector species from infected host species at stage k.

Fig. 1

Flow diagram of Model 1. The dashed, dotted and dash-dotted lines capture the infection routes of vectors from all the infectious hosts at different stages. The planes lines capture the infection routes of hosts from the vectors. We did not display the recovery routes from each I (for k = 1,2…,n − 1 and i = 1,2,…,m) to R.

Table 1

Description of the parameters used in System (2).

Parameters	Description
Λ_h=[Λ_1,Λ_2,…,Λ_m]^T	Vector of recruitment of the hosts;
Λ_v=[Λ_v,1,Λ_v,2,…,Λ_v,p]^T	Vector of recruitment of the vectors;
ai,j	Biting rate of vector j on Host i;
βi,j⋄	Infectiousness of Host i to vectors of species j per biting;
μ_h=[μ_1,μ_2,…,μ_m]^T	Hosts' death rate;
ν_h=[ν_1,ν_2,…,ν_m]^T	Hosts' incubation rate;
α_k=[α_k,1,α_k,2,…,α_k,m]^T	Hosts' total duration at infectious stage k;
η_k=[η_k,1,η_k,2,…,η_k,m]^T	Hosts' recovery rate at stage k;
γ_k=[γ_k,1,γ_k,2,…,γ_k,m]^T	Hosts' progression rate from infectious stage k to k+1;
βk,ji	Vector of species j's infectiousness to Host i at stage k;
μ_v=[μ_1,μ_2,…,μ_p]^T	Vectors' natural mortality rates;
δ_v=[δ_1,δ_2,…,δ_p]^T	Vectors' control-induced mortality rates.

Model (2) describes the dynamics of a Multi-Host and Multi-Vector 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], [26], [33] which consider the dynamics zoonoses and one vector species. Our model extend also the multi-host and multi-vector in [15] which explores and SIR-SI type of host-vector interaction. Particularly, our model consists of 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.

Basic properties and reproduction number

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 where In the following lemma, we derive the basic reproduction number for Model (2) following [16], [43].

The basic reproduction number of Model (2) is given by: where denotes the spectral radius operator.

Following the next generation method [16], [43], the system composed of the infected variables in (2) could be decomposed as the sum of two columns vectors and where represents new infections in each host and arthropod species and that of transitions between classes. By letting and , the Jacobian of and , evaluated at the DFE, we obtain: where and It is worth noticing that the matrix has columns and so, the last row of F has columns, in accordance with the first row of F, as has p columns. The matrix V is Metzler (positive off-diagonal entries) and therefore [4]. Given the particular shape of the matrix V, the matrix-wise entries of are given by the induction relationship: where and . Hence, we can deduce that Hence, the next generation matrix is given by: where and, The basic reproduction number is the spectral radius of the next generation matrix. Hence, by denoting the basic reproduction number of System (2) with n hosts, n infectious stages (in each host) and p vectors, we have: . The matrix is an anti-diagonal block matrix, with and the anti-diagonal entries. Therefore, This proves our claim.  □

The expression of the basic reproduction number can be understood heuristically, using biological interpretations as follows. Indeed, since infection process requires a host species and a vector species, the transmission cycle has two steps: from host species to vector species and from vector species to host species. The first step is captured by while the latter is represented by . More specifically, when the host species are entirely susceptible, that is , for , the transmission rate is . Thus, the transmission rate per infectious vector of all species is . However, the average life-span of infectious vectors, of all species, is units of time. Thus, a single infectious vector, accounting for all species, will produce an average of infected host species. Of these host species, a fraction will survive the latency period to become infectious. The probability of these hosts to survive until a stage l of the infection is: During this time, there is an average of infectious bites by susceptible vector species. Thus, during its infectivity period, an infectious host, of all species, will infect on average: susceptible vectors of different species that bite them. Therefore, the net result of these two steps is, since we are dealing with matrices, the spectral radius of or, equivalently of 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 under certain conditions.

Global stability of equilibria

The disease free equilibrium

The disease-free equilibrium (DFE) always exists and is in Ω. The following theorem gives conditions under which it is globally asymptotically stable (GAS).

The DFE is GAS whenever .

The proof consists of proving that all infected variables of System (1) converge to zero and using the local stability of the DFE when to conclude the GAS of the disease-free steady state. Considering that and we obtain, from System (2), that Therefore, we deduce 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]), where is the stability modulus of . Hence, the trajectories of the auxiliary system whose RHS is that of (3) converge to zero whenever . Since all the variables are positive, by the comparison theorem [40], we conclude that Moreover, this implies that: where ⊘ denotes the Hadamard division. Relations (4) and (5) imply the attractivity of the DFE. Moreover, by [16], [43], the DFE is locally asymptotically stable whenever . We conclude thus that the DFE is GAS whenever .  □

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], [17], [23], [36] to obtain a Lyapunov function in the form of where c is the left eigenvector of the next generation matrix.

Existence and uniqueness of the EE

In this subsection, we explore the asymptotic behavior of System (2) when . 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 [42]. 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.

There exists a unique strongly endemic equilibrium for Model (2) whenever and the Host-Vector connectivity configuration is irreducible.

An endemic equilibrium satisfies the relations: In the following, we express all variables at equilibrium in terms of . From (6), we express all () in terms of , as follows: From (7), it follows that , for any , if and only if . Moreover, we can also express in terms of and since Using the equilibrium relation stemming from the equation of infected vectors (in (6)), the relation , along with (7) and (8), we obtain where and and . Now, we want to express S in terms of . That will make (9) in terms of 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 is an equilibrium of the auxiliary dynamical system: Now, we will show that System (11) has unique equilibrium if the connectivity matrix is irreducible and . To this end, we will use elements from monotone systems and particularly Hirsch's theorem [21]. Indeed, System (11) is monotone if and only if the vector field is, or equivalently whenever the is Metzler. That is, the off-diagonal elements of are positive. We have: since the first term is a diagonal matrix. It follows from (12), that is Metzler as long as is. We have: as . Hence, M and therefore is Metzler. Moreover, since the Host-Vector configuration is irreducible, is irreducible. Thus, the auxiliary dynamical system (11) is strongly monotone. Moreover, it could be seen that the map of is monotonically decreasing. Since and , by Hirsch's theorem ([21], page 55), the trajectories of the auxiliary system (11) either tend to the origin, or else there is a unique equilibrium and all trajectories tend to . The origin is unstable since, for , we have . Indeed, it could be shown easily that . We have shown that if and the host-vector configuration is irreducible, a unique solution of Equation (10) exists. Hence, using (6), (7), we deduce that and , for . We conclude therefore that System (2) has a unique strongly endemic equilibrium whenever and is irreducible.  □

The following subsection consists of investigating the asymptotic properties of this unique strongly endemic equilibrium.

Global stability of the EE

The strongly endemic equilibrium is GAS whenever it exists.

We consider the following Lyapunov function candidate where The coefficients are positive to be determined later. The function 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 along the trajectories of (2) is given by: where and .

Given the expression of , the terms in sum to zero. Similarly, using the expression of , the terms in also sum to zero. Hence, what remains in (14) is: Our first major claim is: Indeed, the equality (★) claims that all linear terms in () in (15) sum to zero. By using the expressions of , we obtain: Substituting this expression of in (★), we obtain: We work the second sum in (★★) out in order to obtain a factor of and therefore compare it with the first one. We have: obtained by replacing the index l by k for convenience. Now, by successively switching the subindices i and r; and using properties of nested sums, we obtain: Thus, showing (★★) is equivalent to show that: This is also equivalent to: Showing (★★★) concludes the proof of (★). We notice that the relation (★★★) is satisfied if are the components of the solution of the linear system , with where and Moreover, we have: Since the Host-Vector connectivity configuration is irreducible, the matrix is also irreducible. Thus, it could be shown easily that . Moreover, where is the cofactor of the i-th diagonal of . Also, for all . Thus, there exists such that . We choose the coefficients of in Lyapunov function are as such.

To show (★★★), we start by its left hand side (LHS): However, implies that, for any i, . Therefore, Thus, using the properties of nested again, becomes which is exactly the right hand side of (★★★).

In summary, Relation (★), obtained through (★★) and (★★★) cancels linear terms in , for , in the expression of given in relation (15). Thus, the latter yields to: By using endemic relations (13) and the expression of the coefficients , we have and Substituting these expressions and switching i by r in the last sum of (17), the expression becomes: Also, notice that: Therefore, Equation (18) implies that: Moreover, since . Thus, Equation (19) implies Finally, by combining the second and third sums in (20), we obtain: The terms and are clearly definite-negative. Note that, taken separately, the two sums and in (21) are not definite negative. Hence, to finish the proof, we will prove that the sum is definite-negative. To do so, we look at the coefficients of the sums from the graph-theoretical standpoint, following the same approach as [17], [18], [19], [23]. Indeed, let be the directed graph that represents the connectivity between the m hosts (including the l stages) and p vectors. Since the Host-Vector connectivity configuration is irreducible, it follows that graph is strongly connected. Recall that, for , are components of the solution of the system , where is given by (16). The matrix is the so-called Laplacian matrix and its associated the graph is strongly connected since is irreducible, since is. Moreover, the solution of is given by the Kirchhoff's matrix tree theorem [7], [32] as follows: where is the set of all spanning trees T of rooted at Host i. Particularly, where is the set of all arcs in T. For our setup, an arc describes an infection arc that starts from Host i, at stage l directed to Host through Vector j. From a modeling standpoint, connects Host i to all vectors (through ) and connect all vectors to all hosts but i (through , with ). By using Cayley's formula [1], [32], each is the sum of terms, each of which is the product of with . These represent the weight of each spanning tree , rooted at Host i.

Now, we investigate what each term of represents in terms of the graph and its spanning trees. Indeed, each term in for all is a weight of a unicyclic graph of a particular length, obtained by adding an arc to a directed tree rooted at Host i, . The obtained unicyclic graph Q has a unique cycle CQ of length . We group all isomorphic cycles as they will have the same coefficients. Hence, we conclude that is the sum over all unicyclic graphs as , where where , (for which ) and represents the cycles that correspond to elements in , and exclusively. Now, each of the three sums in (22) are definite-negative. Indeed, we have: since is a cycle. Also, Similarly, the product of the non constant terms in the last sum (22) is equal to 1 as is a cycle.

Hence, by the arithmetic-geometric mean, is definite-negative for each unicyclic graph Q. It follows from (21), that 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 [17], [18] 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 where is is the Lyapunov function for one host and multiple vectors, is not definite-definite negative for the coefficient .

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 where where And This Lyapunov function has its derivative along the trajectories of the system, definite negative. The corresponding coefficients are determined in the same fashion as the case.

Conclusion and discussions

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], [26], [33] —by incorporating multiple arthropod vector species— and [14], [15] —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 , for m hosts, n hosts' infectious stages and p vectors species. We proved that the disease free equilibrium is globally asymptotically stable whenever (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 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 [15] 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 [17], [18], [23]. 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 consist of constant recruitments.