Modeling ebola virus disease transmissions with reservoir in a complex virus life ecology.

We propose a new deterministic mathematical model for the transmission dynamics of Ebola Virus Disease (EVD) in a complex Ebola virus life ecology. Our model captures as much as possible the features and patterns of the disease evolution as a three cycle transmission process in the two ways below. Firstly it involves the synergy between the epizootic phase (during which the disease circulates periodically amongst non-human primates populations and decimates them), the enzootic phase (during which the disease always remains in fruit bats population) and the epidemic phase (during which the EVD threatens and decimates human populations). Secondly it takes into account the well-known, the probable/suspected and the hypothetical transmission mechanisms (including direct and indirect routes of contamination) between and within the three different types of populations consisting of humans, animals and fruit bats. The reproduction number R0 for the full model with the environmental contamination is derived and the global asymptotic stability of the disease free equilibrium is established when R0andlt;1. It is conjectured that there exists a unique globally asymptotically stable endemic equilibrium for the full model when R0andgt;1. The role of a contaminated environment is assessed by comparing the human infected component for the sub-model without the environment with that of the full model. Similarly, the sub-model without animals on the one hand and the sub-model without bats on the other hand are studied. It is shown that bats influence more the dynamics of EVD than the animals. Global sensitivity analysis shows that the effective contact rate between humans and fruit bats and the mortality rate for bats are the most influential parameters on the latent and infected human individuals. Numerical simulations, apart from supporting the theoretical results and the existence of a unique globally asymptotically stable endemic equilibrium for the full model, suggest further that: (1) fruit bats are more important in the transmission processes and the endemicity level of EVD than animals. This is in line with biological findings which identified bats as reservoir of Ebola viruses; (2) the indirect environmental contamination is detrimental to human beings, while it is almost insignificant for the transmission in bats.


22
T. BERGE, S. BOWONG, J. LUBUMA AND M. L. M. MANYOMBE parameters on the latent and infected human individuals. Numerical simulations, apart from supporting the theoretical results and the existence of a unique globally asymptotically stable endemic equilibrium for the full model, suggest further that: (1) fruit bats are more important in the transmission processes and the endemicity level of EVD than animals. This is in line with biological findings which identified bats as reservoir of Ebola viruses; (2) the indirect environmental contamination is detrimental to human beings, while it is almost insignificant for the transmission in bats.

Introduction.
Since 1976, about 25 outbreaks of EVD have been declared including the latest one in Western Africa (2014)(2015) which was the most devastating one among human populations [18,56]. There are five known Ebola virus strains: ZEBOV, SEBOV, ICEBOV, BDBV, REBOV [18]. Of these, the first four have highly threatened both human and nonhuman primates, causing viral hemorrhagic fever with case fatality rates of up to 90% [18,50]. For instance, the most recent and deadliest Ebola outbreak in Western Africa was caused by the ZEBOV virus strain.
No approved treatments and homologated vaccines are currently available. However, dedicated research efforts have led to the first therapeutic trial with ZMapp [28]. Furthermore, the VSV-ZEBOV vaccine has been found and is still in its third experimental phase [28]. These two remarkable efforts have helped to curtain the recent Western Africa outbreak, even though the latter effort was less decisive (due to its late trial at the end of the West African outbreak) than the former.
Other control efforts (including rehydration, isolation, education of populations at risk, avoidance of consumption of bush meat, practicing of safe burials) have been implemented to stop past Ebola outbreaks. In addition to the threat EVD poses to human health, the negative impact of EVD infection on already threatened animal populations in Africa has come to light and led to a resurgence of efforts to understand the complex life ecology of Ebola virus in nature [26,51]. However, despite considerable efforts, it remains unclear how the EVD is maintained and transmitted in nature, and how the index case (first patient) is infected. Since EVD is a zoonotic-borne disease (transmitted accidentally by direct contact with infected living or dead animals), human epidemics were concomitant with epizootic in great apes [37,38]. Moreover, due to the fact that recent works have provided new evidence that fruit bats might play an important role as a reservoir species of EVD, some intricate and pending questions which are biologically relevant have been raised: 1. Where and how does the index case (first patient) acquire the infection? 2. Do direct transmissions from bats to humans and/or nonhuman primates occur? 3. Which human behaviors expose humans to the risk of contracting EVD from non-human sources? 4. Are fruit bats the only reservoir hosts for Ebola viruses? 5. What are the environmental factors contributing to Ebola virus transmission to human beings and non-human primates from the reservoir species? 6. Can mathematical modeling help to understand and predict the of EVD outbreaks in the future? This paper focuses on the last question. We address it, with the ultimate aim to answer the other five questions. Before moving to the modeling section, there is a need to support and motivate these questions. Since the 1976 Sudan outbreak where there was evidence that the index case was a worker in a cotton factory with evidence of bats at site, many other index cases (from 1994-2001) showed evident contacts with bats and/or consumption of butchered great apes and/or other wildlife meat. The synoptic Table 1 summarizes the sources for contamination of the first patient during EVD outbreaks and highlights his/her contact with bats, dead or butchered wildlife bush meat. The table also talks to the transmission from bat to human, from non-human primate to human, and it illustrates how human behaviors can drive the contamination of EVD. The issue of fruit bats being reservoir for Ebola viruses is supported by the work [37] which has demonstrated that fruit bats bear Ebola viruses and are not affected by the disease. Furthermore, the first four questions are biologically investigated in [26,37,51] where the known, the probable/suspected and the hypothetical direct transmission mechanisms (routes) of EVD are addressed. Regarding the fifth question, the works in [9,10,46,53,60] highlight the indirect environmental contamination route of EVD.
The complexity of the questions raised above is captured in [26,51] for the EVD mechanisms of transmission in a complex Ebola virus life ecology as depicted in Fig. 1 of the disease transmission diagram. This will be reflected in the construction of our model by: (a) Taking into account the well known, the probable/suspected, the hypothetical and the environmental transmission pathways; (b) Involving the interplay between the epizootic phase (during which the disease circulates periodically amongst non-human primate populations and decimates them), the enzootic phase (during which the disease always remains in fruit bat populations) and the epidemic phase (during which the EVD threatens and decimates human populations) of the disease under consideration.
The complexity of the Ebola virus life ecology is clear from the biological studies carried out in [26,51] (see also Ebola virus ecology and transmission therein). From the mathematical point of view, the complexity and challenges we are confronted with and have addressed range from the modeling of the forces of infection, the computation of the reproduction number (see Eq. (8) and Remark1), the computation of the endemic equilibrium, to the use of less standard tools (e.g. decomposition techniques [54]) to investigate the dynamics of EVD.
Note that the recent and former EVD outbreaks have highlighted the importance of human behavior in the transmission process [11,21,41]. For instance, there were evidence of behavioral reaction and self-protection measures: people were scared, they panicked and left care centers, etc... As this important feature calls for a different modeling approach based on "Behavioral Epidemiology" developed in [41], we are busy investigating it in another work, where self-protection measures driven by human behavior is incorporated.
The rest of the paper is organized as follows. In Section 2 we derive the model. The mathematical analysis starts in Section 3 with the basic properties of the model, followed by the computation of the reproduction number R 0 and the establishment of the global asymptotic stability of the disease-free equilibrium. It ends with the investigation of the endemic equilibrium. Section 4 deals with the model without the environmental contamination, while the sensitivity analysis is shown in Section 5. In Section 6, we provide numerical simulations to support the theory and assess the impact of the environmental contamination. Finally Section 8 summarizes our findings and highlights possible extensions. 2. Model formulation. To avoid confusion, bat will not be called "animal". The term animal is reserved for any non-human primate and/or any other wild animal that may be responsible for the transmission of EVD. Note that our model formulation is focused essentially on non-human primates (great apes) as animal species. Motivated by the biological papers [25,26,37,51] regarding the transmission mechanisms and the recent mathematical works [9,53], our model is based on the zoonotic-borne disease setting and takes into account both direct and indirect transmission routes. We distinguish three host populations: humans and animals as end hosts and fruit bats as intermediate reservoir hosts [5,32,37,45,47]. Since the model incorporates the indirect environmental transmission, we add the dynamics of the concentration of free living Ebola viruses in the environment [53,56,60]. The latter compartment should not be considered as an epidemiological class; it is regarded as a pool of Ebola viruses. This pool is supplied by infected humans, infected animals and infected fruit bats with: (i) The presence of carcasses of infected and dead animals in the forest on which some animals can feed [56]. (ii) The manipulation of infected fruit bats and animals hunted by humans for food. (iii) The contaminated fruits harvested for food by humans and primates in the forest [25,26,37]. (iv) The contaminated syringes re-used in health care centers [50,55,58].
(v) The bed linen contaminated by infected human's stool, urine, vomits or sweat in health care centers or in family homes of infected individuals [50,55]. (vi) The bush fruits contaminated by bats droppings [38].
2.1. The variables. As mentioned earlier, the variables include the human, the animal, the bat (reservoir of Ebola viruses) populations and the concentration of free living Ebola viruses in the environment. More precisely, at time t, S h (t), S a (t) and S b (t) denote the susceptible human, animal and bat compartments, respectively. The symbols I h (t), I a (t) and I b (t) denote the infected human, animal and bat compartments, respectively. Since there is an intrinsic incubation period of approximately 2-21 days (on average 8-10 days) of EVD in humans [7,50,55], we introduce E h (t), the exposed (or infected in latent stage) human individuals class. There is a recovered class R h (t) for humans only, because fruit bats are the reservoir of Ebola viruses and nobody cares for infected animals (once infected, they ultimately die). The total human, animal and bat populations at time t are, respectively. The concentration of Ebola viruses in the environment at time t, is denoted by V (t).

2.2.
Main assumptions. The model derivation is based on the following main assumptions: (A 1 ) Infected humans always deposit Ebola viruses in the environment through the routes indicated above, with which susceptible individuals can come into contact. (A 2 ) Human-animal and human-bat contact rates are very low. This assumption is motivated by the fact that the literature does not indicate clearly how the human beings come into contact with either living animals, or the living fruit bats. (A 3 ) Contacts between animals and fruit bats are frequent, occurring during competition for food in the dry seasons where food (specially fruits) is scarce. (A 4 ) Infected fruit bats always contribute to the Ebola virus shedding in the environment. Actually, as reported in [25,26,37,38], the scarcity of food during the dry season, compels bats to deliver allowing their placenta and blood to contaminate fruits and leaves on which some animals (apes, duikers) feed. (A 5 ) The sharing of contaminated food (fruits) by animals (monkeys, duikers) and fruit bats is possible allowing adequate contacts between these two wildlife species. (A 6 ) There are not Ebola-deceased fruit bats since they are a natural reservoir of Ebola viruses.
(A 7 ) Ebola-deceased humans can still infect during unsafe funeral practices, where their corpses are manipulated (e.g. washing, autopsy, dressing-up) as it is the case in many places in Africa [50]. (A 8 ) Ebola-deceased animals can still infect humans, e.g., during manipulation of bats for food [50]. (A 9 ) Clinically Ebola-recovered men (resp. women) still transmit the disease probably through sexual intercourse (resp. breast-feeding) [50]. (A 10 ) Ebola-infected animals do not recover since nobody cares about them.
2.3. The equations. The works in [26,51] amongst others on the complex Ebola virus life ecology lead to the dynamics flowchart in Fig. 1 which in turn gives the following system of nonlinear ordinary differential equations: The force of infection acting on humans is, (2) It is the sum of the contributions below: • The human-to-human force of infection which gathers the three contamination processes by: • The animal-to-human force of infection • The bat-to-human force of infection

27
• The environment-to-human force of infection Similarly, the force of infection within the animal population involves the following three contributions: • The animal-to-animal force of infection • The bat-to-animal force of infection • The environment-to-animal force of infection Finally, the force of infection in bat's population is modeled by and consists of two contributions from: -Within bat adequate contacts -The contact with the environment Note that we have considered the infection through contact with environmental free Ebola viruses. As it is the case for most models involving free-living pathogens in the environment [3,4,8,12,17,49,53,59], the environmental-related forces of infection, λ hv , λ av and λ bv are modeled using Michealis-Menten or Holling type II functional responses. The constant K represents the minimum amount of viruses in the environment capable of ensuring 50% chance of contracting the disease. Moreover, for the sake of simplicity, the contribution of dead animals and humans to environment is modeled by relating them to infective individuals. This simplification avoids the need to add compartments counting number of deceased humans and animal carcasses in the environment.
Finally, the last equation of system (1) describes the dynamics of Ebola viruses with shedding from humans, animals and bats. The parameters used for system (1) and their biological interpretations are giving in Table 2. For notational simplifications let, , and Φ a = 1 + ξ a ν a δ a .
With this notation it is easy to check that system (1) interconnects the following sub-models: • The human population sub-model • The animal population sub-model • The bat population sub-model It should be emphasized that not much is known about the bat-to-bat EVD transmission. However, knowing that fruit bats settle or congregate for rest or sleep (they live in colony), it is acceptable to assume that direct bat-to-bat contact is the main route of transmission and can be modeled by mass action incidence.
• The evolution of free-living Ebola viruses in the environment which is modeled by Eq. (1e). The model presented her is a new in many respects. It extends the existing works [1,9,16,24,35,36,44,53] in the sense that it considers the three phases of disease: the epizootic cycle in animals, the enzootic phase in fruit bats and the epidemic phase in humans. In particular, the novelty of our model is clear from the most recent work [9], where direct human-to-human transmission was considered and all other sources (e.g. consumption of bush meat, manipulation of fruit bats, indirect environmental contamination) were encompassed in a constant recruitment of Ebola viruses in the environment. The model developed here enriches the latter by modeling this recruitment through consideration of the complex Ebola virus life ecology, where animals and bats are explicitly involved in the EVD transmission cycle. It is therefore understandable why throughout this paper, we refer to system (1) as the full model. 3. Theoretical analysis of the full model.

Positivity and boundedness of solutions.
For the EVD transmission model (1) to be epidemiological meaningful, it is important to prove that all state variables are non-negative at all time. That is, solutions of the system (1) with non-negative initial data will remain non-negative for all time t > 0. (1) are non-negative for all t > 0, when it exists.
, then for any initial condition such that where ρ(t) = exp t 0 (λ h (p) + µ h )dp > 0 is the integrating factor. Hence, integrating this last relation with respect to t, we have so that the division of both side by ρ(t) yields The same arguments can be used to prove that S a (t) > 0,  Furthermore

Bats (Reservoir) Animals
Finally, using the fact that This completes the proof.
Combining Theorem 3.1 with the trivial existence and uniqueness of a local solution for the model (1), we have established the following theorem which ensures the mathematical and biological well-posedness of system (1) (see [9], Theorem 3.3).
To compute the basic reproduction number of the model, we use the standard method of the next generation matrix developed in [2,8,19,20]. We separate the infected states (E h , I h , R h , I a , I b , V ) form the uninfected states (S h , S a , S b , ). Let F and W be the vectors representing the new and transported cases into the infected states, respectively. Thus The Jacobian matrices F of F and W of W evaluated at the DFE are W is a lower triangular and invertible matrix. Thus, thanks to [8], R 0 is obtained as the maximum eigenvalue of the positive matrix F W −1 , where , , .
Since zero is an eigenvalue for F W −1 of multiplicity 3, simple algebraic matrix properties show that its non vanishing eigenvalues are those of the (3 × 3) matrix Therefore, where for a square matrix M , ρ(M ) denotes its the spectral radius. Based on some of the realistic assumptions stated in subsection 2.2, the remark below gives the explicit formula of the basic reproduction number R 0 in special cases.
Remark 1. In some cases, the explicit formula for R 0 are straightforward as thrived below.
1. In the sub-model with only human population dynamics and environmental transmission [53], the basic reproduction number denoted by R 0,hv is If in addition, the indirect transmission is neglected (i.e. β hv = 0), then the basic reproduction number reduces to We emphasize that the basic reproduction number giving by (10) is suitable for comparison with the basic reproduction numbers for some existing EDV models. For instance, looking at the expression Φ h , (10) is the sum of three contributions: (i) The contribution from infected human individuals, i.e. [53] for more details in such comparisons). (iii) The contribution from Ebola-deceased individuals, i.e. β hh ξ h ν h γf ω (µ h +ω)(µ h +γ) . With the above decomposition of R 0,h , it is straightforward that the environmental contamination increases its range, such that R 0,h is always larger than the basic reproduction numbers for the few existing SEIR classical models with standard incidence forces of infection.

2.
In the sub-model without animal's population dynamics, the basic reproduction number denoted by R 0,hbv is 3. In the sub-model without bat's population dynamics, the basic reproduction number denoted by R 0,hav is 4. Suppose β bv = 0, then R hbv 0 = R abv 0 = 0 and the basic reproduction number is 5. If β bv = β av = 0, then R 0 simply becomes where R 0,hv is given in (9) and are the intra-specific basic reproduction numbers of the animal's population and bat's population without the environmental transmission, respectively and R 0,hv is the intra-specific basic reproduction number in the human population with the environmental transmission.
The epidemiological implication of Lemma 3.3 is that EVD can be eliminated from the community when R 0 < 1 and the initial sizes of the sub-populations in the model are in the basin of attraction of the DFE P 0 . But, for the disease to be eliminated independently of the initial sizes of sub-populations, the global asymptotic (GAS) stability of the DFE must be established when R 0 < 1. This is the substance of the following theorem.
Proof. Let x = (E h , I h , R h , I a , I b , V ) and y = (S h , S a , S b ) be the disease compartments (infected) and uninfected states, respectively. Then system (1) can be re-written in the form where F and W are giving above, We then consider the following auxiliary linear subsystem from (16): From Theorem 2 in [20], we have R 0 < 1 ⇐⇒ σ(F − W ) < 0, where, for a square matrix M , σ(M ) denotes its stability modulus. So, when R 0 < 1, the eigenvalues of F −W all have negative real parts. Therefore, non-negative solutions of (17) are such that lim t→+∞ x = 0, or equivalently lim t→+∞ E h = lim t→+∞ I h = lim t→+∞ I a = lim t→+∞ I b = lim t→+∞ V = 0. By the standard comparison principle [33,48] and the non-negativity of x, non-negative solutions of (1) satisfy lim t→+∞ E h = lim t→+∞ I h = lim t→+∞ I a = lim t→+∞ I b = lim t→+∞ V = 0. Therefore, since lim t→+∞ x = 0, system (1) is an asymptotically autonomous system [14] (Theorem 2.5) with the limit system: It is straightforward that the linear system (18) has a unique equilibrium given by (S 0 h , S 0 a , S 0 b ) which is globally asymptotically stable. This completes the proof. Remark 2. To extend this global result in the case when R 0 = 1 (which we do not address here), the construction of a suitable Lyapunov function and the use of LaSalle' Invariance Principle are necessary.

3.3.
Existence of endemic equilibrium of model (1). In this section, we investigate the existence of equilibrium points other than the disease free equilibrium, namely possible boundary equilibrium points and interior equilibria. First of all, let us give some useful remarks.
Assume that an equilibrium is such that I b = 0, then from Eqs. (1h) and (1i), V = 0, and I a = 0. Using these statements in Eqs. (1b), (1c) and (1d) gives Similarly, if an equilibrium of (1) is such that I a = 0, then from Eq. (1e), I a = V = 0 and from Eq. (1i), one has I h = 0. Replacing these values in Eqs. (1b), (1c) and (1d) yields I h = E h = R h = 0 and the said equilibrium point is disease free as well.
Obviously, if an equilibrium of (1) is such that V = 0, then from Eq. (1i), I h = I a = I b = 0. Introducing these in Eqs.(1b), (1c) and (1d) leads to the E h = R h = 0, and once more, the corresponding equilibrium is disease free. All in all, the only boundary equilibrium point for system (1) where the disease absent in the human population is the disease free equilibrium.
Conversely, assume the human population is disease free, then the free virus concentration V = 0, and from Eq.(1i), we have I a = I b = 0. Thus, the full system is disease free. Note that, the non existence of boundary equilibria is due to the fact the disease transmission "one way" (that is, from animals and bats to humans and not the other way round).
As a consequence, we have proven the following result: Lemma 3.5. System (1) has no other boundary equilibrium than the disease-free equilibrium.
This lemma is very important as it excludes the possibility for the full model (1) to exhibit non trivial boundary equilibrium points. This suggests that the full model could have exactly one interior (endemic) equilibrium with the disease being present in all the populations under consideration. This, together with the existence and uniqueness of interior equilibrium for some system (1)-related sub-models [9], motivates the following conjecture that we make. Conjecture 1. Assume that R 0 > 1 for system (1). Then there exists a unique interior (endemic) equilibrium.
The stability of the endemic equilibrium will be shown numerically at a later stage.
Remark 3. Actually, Conjecture 1 could be addressed in a separate work, by reducing the finding of equilibria to a fixed-point problem and apply a suitable fixedpoint theorem for a multi-variable and sub-linear function for monotone dynamical systems or for systems of ordinary differential equations which generate an order preserving flow [29,30,48].
In order to investigate the effects of the environmental contamination on the transmission of EVD, it is reasonable to consider the sub-model of system (1) without the environment compartment. Note that, even though the full model (which couples many subsystems) could have two equilibria (namely, the disease free and the interior equilibria), it is not obvious that all its sub-models will exhibit the same property since the coupling can reduce or increase for example the number of equilibrium points. Below is the sub-model involving only human, animal and bat populations, but excluding the environment influence.
of infection. The model in this setting reads therefore: The corresponding basic reproduction for this model is easily computed as where, Actually, R 0,h , R 0,a and R 0,b are the intra-specific basic reproduction numbers for human, animal and bat sub-populations given earlier by Eq. (10) and Eq. (15), respectively.
Remark 4. Note that the non-negative matrix G in (8) has diagonal entries R hhv 0 , R aav 0 and R bbv 0 that are larger than R 0,h , R 0,a and R 0,b , respectively. This, together with the fact that the spectral radius of a non-negative matrix is an increasing function of its entries yields Consequently, the indirect environmental contamination enhances the transmissibility of EVD, and thus it increases the epidemic/endemic level of the disease.
The dynamics of sub-system (19) are confined in the biological feasible compact set subset of Γ hab given by Sub-model (19) is triangular. Indeed the variables (S h , E h , I h , R h ) do not appear in the last fourth equations. Moreover, the variables (S a , I a ) do not appear in the last two equations. Therefore, in order to address the long run behavior of system (19), we shall make use of the decomposition techniques by applying repeatedly the following theorem [54].
1. If x * is globally asymptotically stable (GAS) in R n for subsystem and if y * is GAS in R m for the subsystem dy dt = g(x * , y), then (x * , y * ) is (locally) asymptotically stable for system (23). 2. Moreover, if all the trajectories of system (23) are forward bounded, then (x * , y * ) is GAS for the system (23).
In order to apply Theorem 4.1, we need to study some sub-systems of Eq. (19).

4.1.
Dynamics of the bat sub-model (19g)-(19h). We give here the long run behavior of the sub-model by establishing the global stability of its equilibrium points. Obviously this system (24) has two equilibria; namely, the disease-free equilibrium P 0 b = (S 0 b , I 0 b ), whose coordinates are and the endemic equilibrium which exists whenever R 0,b > 1 and whose components are Obviously the dynamics of sub-model (24) are confined in the biological feasible compact set Looking at model (24) in which the mass action law is applied, it is standard to deduce from the Lyapunov-LaSalle techniques its asymptotic behavior summarized in the result below.
Proposition 1. The following statements hold: • If R 0,b ≤ 1, then the disease-free equilibrium P 0 b for subsystem (24) is GAS. It is unstable whenever R 0,b > 1.
• If R 0,b > 1, then the equilibrium P b is GAS.

4.2.
Dynamics of the animal sub-model (19e)-(19f) with the bat population at equilibrium points. Here, we consider the subsystem (19e)-(19f) The dynamics of sub-model (26) are confined in the biological feasible compact set Ω a = (S a , I a ) ∈ R 2 + : S a + I a ≤ Λ a µ a .
Using the global asymptotic stability of equilibria for subsystem (24), the following two subsystems will be considered.

4.2.1.
Dynamics of (26) with the bat population at equilibrium point P 0 b . The variables S b and I b are substituted in (26) by their corresponding values at the diseasefree equilibrium P 0 b . This leads us to the system Direct calculations show that model (27) has two possible non-negative equilibrium states: the disease-free equilibrium P 0 a = S 0 a = Λ a µ a , I 0 a = 0 and a unique endemic equilibrium P a = S a , I a , with .
We have the following proposition: Proposition 2. The following statements are satisfied.
• If R 0,a ≤ 1, then for the subsystem (27), the disease free equilibrium P 0 a is GAS. It is unstable whenever R 0,a > 1.
• When R 0,a > 1, there exists an unique endemic equilibrium P a which is GAS.
Proof. The GAS of the disease-free equilibrium P 0 a is established using the quadratic Lyapunov function V 0 (S a , I a ) = 1 2 I 2 a . The directional derivative of V 0 towards the vector field given in the right-hand side of (27) isV Thus,V 0 ≤ 0 whenever R 0,a ≤ 1. Observe thatV 0 (S a , I a ) = 0 if and only if, either I a = 0 or R 0,a = 1 and I a = 0. In both cases, it is easy to see that the largest invariant set contained in V 0 (S a , I a ) = 0 is reduced to the disease-free equilibrium P 0 a . Hence, by LaSalle Invariance Principle [34], P 0 a is GAS.
For the GAS of the endemic equilibrium P a , we propose the following Lyapunov function candidate.
Q(S a , I a ) = (S a + I a ) − (S a + I a ) − (S a + I a ) ln S a + I a S a + I a defined in the set {(S a , I a ) ∈ Ω a : S a > 0, I a > 0}. The positive constant k will be determined shortly. Since (S a , I a ) is an equilibrium of (27) we have With this in mind and the fact that the directional derivativeQ of Q towards the vector field given in the right-hand side of (27) iṡ Choose the constant k such that 2µ a + δ a − kβ aa Φ a I a S a + I a = 0, or equivalently k = (2µ a + δ a )S a + I a β aa Φ a I a .
Thus, the directional derivative of Q becomeṡ from which we can see clearly thatQ < 0 except at the endemic equilibrium where it is zero. Therefore, P a is GAS.

4.2.2.
Dynamics of (26) with the bat population at equilibrium point P b . Submodel (26) is considered when the bat population is at endemic state P b . That is the variables S b and I b are replaced in (26) by S b and I b , respectively. This gives to the following subsystem: Since I b > 0, (28) has no disease free equilibrium. We study the existence of endemic equilibria. Let D b be defined by E a = S a , I a is an equilibrium point of (28) if and only if Set Then, from (30) and (31), we have Substituting (32) into (31) yields From this latter expression, we derive the quadratic equation (33) Denote the discriminant of Eq. (33) by Then, the unique positive root of Eq. (33) is Thus, the components of the unique equilibrium point E a , obtained by substituting (35) into (32) are Proposition 3. The endemic equilibrium E a for subsystem (28) is GAS whenever R 0,a ≤ 1.
Proof. We first establish the LAS of E a when R 0,a ≤ 1. Let J( E a ) be the Jacobian matrix at any equilibrium point E a of (28). We have Since E a is an equilibrium point, it can be shown that the trace of J( E a ) is which is positive whenever R 0,a ≤ 1. This proves the LAS of E a .
Secondly, we prove the global attractiveness of E a . To achieve this, we use the Dulac criterion [27] to rule out the existence of periodic solutions. Consider the Dulac function g(S a , I a ) = 1 I a , defined on the connected region 0, Λ a µ a × 0, Λ a µ a containing the interior of Ω a .
Let X(S a , I a ) = (X 1 (S a , I a ), X 2 (S a , I a )) T be the right hand side of (28). It is easily shown that Hence, by Dulac's criterion, there is no periodic solution in the interior of Ω a . Hence E a is GAS whenever R 0,a ≤ 1. This completes the proof.
Model (40) is well posed mathematically and biologically in the compact set System (40) has two equilibrium points, namely the disease-free which always exists and the endemic equilibrium point E h = S h , E h , I h R h , which exists whenever the human intra-specific basic reproduction number Then the components of E h are: .

MODELING EBOLA VIRUS DISEASE IN A COMPLEX VIRUS LIFE ECOLOGY 43
The asymptotic behavior of model (40) is completely described by The following statements hold true: (1) The disease-free equilibrium point P 0 h is GAS when R 0,h ≤ 1 and unstable if R 0h > 1.
(2) The endemic equilibrium point E h is GAS whenever R 0h > 1.
Proof. The first item is established using the Lyapunov function The Lie derivative of L h with respect to the vector field given by the right hand side of (40) iṡ Indeed, suppose I h = 0, then replace it in the first, second and fourth equations of (40) and solve. One h . The conclusion for the GAS of P 0 h follows by LaSalle's Invariance Principle.
As for the second item of Theorem 4.2, the proof of the global asymptotic stability is quiet long and challenging. We refer the interested reader to [39,61] where the proof is provided using a geometrical approach [40].
Thanks to Theorem 4.1 and combining Proposition 1, Proposition 2 and Proposition 4.2 we are now able state the first main result regarding the asymptotic behavior of the environmental-free model (19). Theorem 4.3. For system (19), the following statements hold true: 1. If R 0,h ≤ 1, R 0,a ≤ 1, and R 0,b ≤ 1, then P 0 h , P 0 a , P 0 b is GAS . 2. If R 0,h > 1, R 0,a ≤ 1, and R 0,b ≤ 1, then P h , P 0 a , P 0 b is GAS.

4.3.2.
Dynamics of (39) with the bats and animals at steady state P a , P 0 b . The model under consideration in this section is where the constant D a is giving by This model is obviously a dynamical system on the biological feasible domain given in (41). Due to the fact I a > 0, model (42) cannot exhibit a disease free equilibrium. We shall therefore focus on the existence and stability of possible endemic equilibrium points. Let E * * h = (S * * h , E * * h , I * * h , R * * h ) be an equilibrium of (42). Set From (45), we have , .
(46) Putting the above expressions (46) in (44), yields the following quadratic equation where the constant coefficients a h , b h , c h are: Since, a h > 0 and c h < 0, the quadratic equation (47) has one positive solution λ * * h to which corresponds the unique endemic equilibrium point E * * h of subsystem (42) giving by (46). Similar arguments to those in [39,61] can be used to prove the global asymptotic stability of the unique endemic equilibrium point E * * h of subsystem (42). Hence, thanks to Theorem 4.1, we can give the second main theorem of this section. for subsystem (42) is GAS.

4.3.3.
Dynamics of (39) with the bats and animals at steady state E a , P b . Replacing I a , I b in subsystem (39) above, yields the following system where D a is the constant defined by The theoretical analysis of subsystem (49) is similar to that of subsystem (42).
We summarize the existence of the four equilibria of system (19) and their stability properties in the following table.

Equilibria
Conditions of existence Stability Table 3. Existence, conditions for existence and stability of equilibria.
5. Sensitivity analysis. We carried out a sensitivity analysis. This allows to identify the parameters that are most influential in determining population dynamics [42,43]. A Latin Hypercube Sampling (LHS) scheme [15,43] samples 1000 values for each input parameter using a uniform distribution over the range of biologically realistic values, listed in Table 7. Using system (1), 1000 model simulations are performed by randomly pairing sampled values for all LHS parameters. Outcome measures are calculated for each run: Exposed individuals to EVD, infected individuals, virus concentration, infected animals and infected bats. Partial Rank Correlation Coefficients (PRCC) and corresponding p-values are computed. An output is assumed sensitive to an input if the corresponding PRCC is less than −0.50 or greater than +0.50, and the corresponding p-values is less than 5%. The results are displayed in Table 4, Table 5 and Table 6. From these tables, it can be seen that the effective contact rate between humans and fruit bats and the bat mortality rate are the most influential parameters on the latent and infected human individuals.
Parameters  Table 4. PRCCs of full model's parameters 6. Numerical simulations. In this section, we give numerical simulations that support the theory presented in the previous sections. The simulations are produced by MatLab. While the parameters values for human-to-human transmission are mostly taken from [50,57], we have proposed almost all the parameter values regarding animal-to-human, bat-to-human, environment-to-human, environment-toanimal, environment-to-bat, animal-to-animal, animal-to-bat and bat-to-bat transmission mechanisms.
6.1. General dynamics. We numerically illustrate the asymptotic behavior of the full model and the sub-model without the environmental contamination. The GAS of the disease-free equilibrium P 0 demonstrated in Theorem 3.4 and the existence and stability of a unique endemic equilibrium as stated in Conjecture 1 for the model with the environmental contamination are numerically shown on Figure 2 and Figure 3, respectively. However, Figure 4 further suggests the GAS of (P h , P 0 a , P 0 b ).   Table 6. PRCCs of model's parameters without animals Figure 5 illustrates the GAS of (P 0 h , P 0 a , P 0 b ) for the free-environmental contamination sub-model (19) as established in Theorem 4.3. Figure 6 supports the stability of (E * * h , P a , P 0 b ) as shown in Theorem 4.4, while Figure 7 illustrates the stability of (E * * * h , E a , P b ) as shown in Theorem 4.5.  Table 7. Baseline numerical values for the parameters of system (1) 6.2. Impact of the contaminated environment on the infected level of EVD. We numerically assess the impact of the contaminated environment on the severity/endemicity of EVD, as well as the effects of bats and animals on the long run and severity/endemicity of EVD. Figure 8 and Figure 9 show the increasing behavior of the full model-related infected component with respect to the indirect effective contact rates β hv , β av and β bv . This highlights the detrimental role of the contaminated environment on the transmission dynamics of EVD. Moreover, one observes from Figure 8 that, the infected bats are not influenced by the environmental contamination, which suggests that the indirect route of transmission for bats can be neglected. Similarly Figure 9 shows that the infected humans and bats are not significantly influenced by the incorporation of animals species in the model development. This suggests and confirms in some sense (see [37]) the fact that the bats are reservoir of Ebola viruses. Therefore, it may be more suitable to build a model involving only two hosts: namely, humans and fruit bats. To illustrate the effects of bats and animals on the dynamics of EVD, Figure  9 is the simulation of the model (1) with and without bats on the one hand and with and without animals on the other hand. It can be seen that the incorporation of bats increases significantly the endemic level of EVD in the human population, while the involvement of animals does not.
7. Conclusion. The main purpose of this paper was to build and analyze a mathematical model for the transmission dynamics of EVD in a complex Ebola virus life ecology. We have then developed and analyzed both theoretically and numerically a new model by taking into account the known, the probable/suspected and the hypothetical mechanisms of transmission of EVD [26,37,51]. The proposed model captures as much as possible the essential patterns of the disease evolution as a three cycle transmission process in the following two ways: 1. It involves the interplay between the epizootic phase (during which the disease circulates periodically amongst non-human primates populations and decimates them), the enzootic phase (during which the disease always remains in fruit bats population) and the epidemic phase (during which the EVD threatens and decimates the human beings). 2. It includes the direct transmission mechanism between and within the three different types of populations which are humans, animals and fruit bats, as well as the indirect route of infection through a contaminated environment.  Figure 7. Stability of (E * * * h , E a , P b ) when Λ h = 10, Λ a = 10, Λ b = 10, µ h = 0.033, µ a = 0.2, µ b = 0.29, δ a = 0.05, α h = α a = α b = 0.95, f = 0.50, β hh = 0.3, β bb = 0.05, β hb = β ha = 10 −8 , β ab = 0.05, β aa = 0.2 (so that R 0,h = 1.7269, R 0,a = 0.8074, R 0,b = 3.1250).
More precisely, we have extended and enriched the few existing SEIR-type human models for EVD with five additional compartments which model the direct transmissions within/between animal and fruit bat populations as well as the environmental indirect contamination. In this double setting of direct and indirect transmissions, our major findings from the theoretical, numerical and computational point of view read as follows: From the theoretical perspective, our results are two-fold: • For the full model with the environmental contamination, we have computed the basic reproduction number R 0 and used it to prove the global asymptotic stability of the disease free equilibrium, whenever R 0 < 1. Furthermore, when R 0 > 1 the existence and the stability of an endemic equilibrium is investigated and conjectured.
• The sub-model without the environmental contamination exhibits one globally asymptotically stable disease-free equilibrium whenever the host-specific basic reproduction numbers, R 0,h , R 0,a and R 0,b are less than or equal to the unity. We have also shown that a such sub-model has three additional endemic equilibria which are all globally asymptotical stable.
From the numerical and computational point of view, the following three facts were addressed:   • In order to assess the role of a contaminated environment on the spreading of EVD, the infected human component resulting from sub-model without the environmental contamination was compared with that of the full model. Similarly, we have considered the sub model without animals on the one hand, and the sub model without bats on the other hand, and found that bats influence more the dynamics of EVD than animals. This is probably because almost all EVD outbreaks were due to consumption and manipulation of fruits bats.
• Global sensitivity analyses were performed to identify the most influential model parameters on the model variables. It shows that the effective contact rate between humans and fruit bats and the bat mortality rate were the most influential parameters on the latent and infected human individuals. This is probably because almost all EVD outbreaks were due to consumption and manipulation of fruits bats.
• Numerical simulations, apart from supporting the theoretical results and the existence of a unique global stable endemic equilibrium for the full model (when R 0 > 1), have further suggested the following two important statements: (1)-fruit bats are more important in the transmission processes and the endemicity of EVD than the animal species. This is in line with biological findings through which fruit bats were identified as the reservoir of Ebola viruses. (2)-the indirect environmental contamination is detrimental to human beings and is almost insignificant for the transmission in bats. From all these investigations, we believe that a more realistic mathematical model will involve only human beings and fruit bats.
Despite the high level of generalization and complexity of our work, it still offers many opportunities for extensions. Theses include: (i) The incorporation of the Ebola-deceased compartments to better capture the transmission mechanisms of EVD during funerals [9]. (ii) The incorporation of the transmission in health care centers in which medical staff can be infected as well [58]. (iii) The incorporation of patches to account for the internationalization of EVD as it is the case in Western Africa [13,23,31]. (iv) The modeling of multi-species transmission mechanism in the case where the same region is threaten by more than one Ebola virus strain. (v) The incorporation of the human behavior. For instance, there were evidence of behavioral reaction and self-protection measures: people were scared, they panicked and left care centers, etc...Thus the need to fill the gap of lack of modeling human behavior [11,21]. This important feature calls for a modeling approach based on "Behavioral Epidemiology" developed in [41], which we are already addressing in another work which takes into account self-protection measures driven by human behavior. (vi) The modeling of some optimal control strategies such as vaccination, isolation, quarantine, treatment, early detection, environmental decontamination.