A simple mathematical model for Ebola in Africa

ABSTRACT We deal with the following question: Can the consumption of contaminated bush meat, the funeral practices and the environmental contamination explain the recurrence and persistence of Ebola virus disease outbreaks in Africa? We develop an SIR-type model which, incorporates both the direct and indirect transmissions in such a manner that there is a provision of Ebola viruses. We prove that the full model has one (endemic) equilibrium which is locally asymptotically stable whereas, it is globally asymptotically stable in the absence of the Ebola virus shedding in the environment. For the sub-model without the provision of Ebola viruses, the disease dies out or stabilizes globally at an endemic equilibrium. At the endemic level, the number of infectious is larger for the full model than for the sub-model without provision of Ebola viruses. We design a nonstandard finite difference scheme, which preserves the dynamics of the model. Numerical simulations are provided.


Introduction
In 2014, an outbreak of Ebola virus (Ebola) decimated many people in Western Africa. With more than 16,000 clinically confirmed cases and approximately 70% mortality cases, this was the more deadly outbreak compared to 20 Ebola threats that occurred since 1976. In almost all the outbreaks, the index case (first patient) became infected through contact with infected animals (hunted for food), such as fruit bats and primates (ape, monkey, and chimpanzee). This suggests that the virus can be spread by indirect contact [13]. As reported in [7], a non-negligible percentage of the Ebola-Zaire virus type survived after 14 days at 4 • C on glass and plastic (10%) and on surfaces (3%). Moreover, 0.1-1% of Ebola virus particles remained active for up to 50 days at 4 • C [31]. This survival of the virus in the environment, due to poor hygienic and sanitary conditions, is probably another source of Ebola infection in many places in Africa. In Africa, and particularly in the regions that were affected by Ebola outbreaks, people live close to the rain-forests, hunt bats and monkeys and harvest forest fruits for food [23,24]. As part of their tradition and customs, Africans are warmhearted to the extent that even a contagious disease would not stop them from caring for their relatives at home, kissing themselves and shaking hands. Furthermore, during funerals, they wash and dress up their deceased relatives. They share, without proper washing, clothes of their deceased relatives. The huge gatherings of people from surrounding villages, towns and cities are most likely among the major factors for the quick spread of the Ebola virus infection in the region.
The above-mentioned considerations raise the following research question: Can the consumption of contaminated bush meat, the funeral practices, and the environmental contamination explain the recurrence and persistence of Ebolavirus disease (EVD) outbreaks in Africa? Naturally, this question should be coupled with the well-known direct transmission route, which involves contact with: (1) blood or body fluid (including but not limited to urine, saliva, sweat, feces, vomit, breast milk, and semen); (2) objects (e.g. clothes, bedding) that have been contaminated with body fluid. To address this question, we develop a simple deterministic model that captures the African practices. The model incorporates both the direct (human-to-human) and the indirect environmental (environment-to-human-to-environment) transmission routes [17,23,24,38]. It is an extended susceptible-infectious-recovered-deceased (SIRD) model with an additional environmental compartment (P) referred to as the pool of Ebola viruses.
The few existing works [1,2,10,15,16,19,21,22,30,32,35,37] on the mathematical modelling for the transmission of the Ebola virus focus only on the human population and direct transmission. As a matter of fact, almost all the mathematical studies for the transmission dynamics of EVD outbreaks are described in classical settings, such as SI model [37], SIR model [16,32], SEIR model [2,10,22], SEIRD model [15,30], or SEIRHD model [1,19,21] where (S) stands for Susceptible, (E) for Infected, (I) for Infectious, (H) for Hospitalized, (R) for Recovered/Removed, and (D) for Dead/Deceased. In the Central Africa case, these studies [10,21,22,30] aimed mostly at estimating the epidemiological parameters and approximating the basic reproduction number. Thus, the earliest work [10] published in 2004, used SEIR model to estimate R 0 at 1.83 for Congo and 1.34 for Uganda. Later in 2006, an SEIR model [22] led to estimates R 0 1.33 for Congo and R 0 1. 35 for Uganda. A different work based on an SEIRD model [21] and conducted in 2007 gave the approximation R 0 2.7 for both the 1995 Congo and the 2000 Uganda outbreaks. The most recent work [21] for the Central Africa outbreaks was done in 2013. There, the authors used a Markov Chain Monte Carlos method to estimate the parameters of their SEIRD model and approximated the pre-intervention R 0 at 2.1 and the post-intervention R 0 at 0.1339 for the 1995 Congo outbreak.
The emergence from 2014 of yet another Ebola outbreak in Africa, namely in Western Africa has created a surge in the mathematical modelling of the EVD. A brief comparison of the modelling approaches used with those of the Central Africa outbreaks is in order. We restrict ourselves to three of the main features of the disease outbreak in Western Africa.
Let us talk first about the demographic dynamics. At the beginning of the Western Africa outbreak, nobody expected the threat to last for more than two years. This constitutes a sufficiently long period of time for the demographic process to be allow to occur. However, none of the mathematical models investigated in 2014-2015 included vital dynamics [2,10,15,16,35]. It is only towards the end of 2015 that modellers, realizing that the outbreak was predicted to persist for several months should nothing be done to stop it  [35], started incorporating the demographic dynamics in their works [1,19]. Secondly, the Western Africa outbreak also distinguished itself by its internationalization which resulted in the development of meta-population models [14,19]. Similarly to the works proposed in the Central Africa case, there has been a great deal on the estimation of the basic reproduction number R 0 for most of the models studied from 2014, as summarized in Table 1.
Since the first outbreak in 1995 in Congo, it was established that the index case (firstinfected individual), either ate contaminated bush meat hunted for food or had contact with fruit bats [9,23,24]. Equally, in the Western Africa outbreak, it was identified that the first case had contact with fruit bats and that the disease started in a small village near rainforest [9,32,38]. More importantly, the survival and the persistence of the Ebola viruses in the environment as well as the environmental transmission of EVD have been recently demonstrated [7,9,31,32,38]. Thus, the environment on the transmission dynamics of EVD is the third main feature of the disease. This important characteristic has to the authors' knowledge not been captured in the modelling of EVD.
The main purpose of this paper is to incorporate (in a simple manner) the indirect environmental transmission on the dynamics of EVD and to assess the effect of such a feature on the long run of the disease.
Our model extends and enriches several of the existing works in four main directions as follows: (i) We incorporate the transmission of deceased individuals during funerals. A compartment for dead individuals is explicitly incorporated in [1,15,19,22,30]. (ii) We include the infection through the contaminated environment resulting from African practices, hospitality and poor hygienic conditions. The consideration of this important feature is new. (iii) Unlike all the existing models in the literature, we include the provision or recruitment source of Ebola virus linked to the consumption of bats, hunted meat and fruits from rain-forests. (iv) We allow the demographic process (vital dynamics) to take place during the EVD outbreaks. This aspect has only been considered recently [1,19].
The model is carefully analysed theoretically and numerically. From the theoretical point of view, we show in the following precise manner that the severity of the disease increases with the recruitment/provision of Ebola viruses and the disease dies out in the absence of such recruitment/provision as well as in the absence of shedding and manipulation of deceased individuals: (1) The full model has only one endemic equilibrium, which is LAS while in the absence of shedding or manipulation of deceased individuals, this equilibrium is globally asymptotically stable (GAS). (2) For the model without recruitment/provision of Ebola viruses, the disease-free equilibrium is GAS when the threshold quantity R 0 is less than or equal to unity. When the said threshold is greater than one, the endemic equilibrium is LAS in the generic case and GAS in the absence of shedding and manipulation of deceased individuals. (3) At the endemic level, for both cases in items (1) and (2) above, the number of infected individuals when there is provision of Ebola viruses is larger than the corresponding number of infected individuals in the absence of such provision.
As it is the case for most systems of differential equations that govern real-life situations, the deterministic model proposed in this work can unfortunately not be solved explicitly by analytic techniques. It is therefore essential and vital to design reliable numerical methods that capture the essential properties of the continuous model. It is well known that classical methods such as the Euler and the Runge-Kutta methods do not serve the purpose [3,4,12]. Thus, from the numerical point of view, we apply the nonstandard finite difference (NSFD) approach initiated by Mickens three decades ago, given its potential to replicate the dynamics of continuous models in many applied areas such as ecology and epidemiology [4,11,18,[26][27][28][29]. In this new discrete framework, we conduct numerical tests based on recent works reporting on Western Africa outbreaks [2,15,35]. Our NSFD scheme, designed on the basis of an innovative use of Mickens rules on the nonlocal approximation of nonlinear terms and complex denominator function of discrete derivatives, preserves the positivity, the boundedness of solutions as well as the number and the stability of equilibrium points. To the authors' best knowledge, the NSFD approach has never been used for the discrete models for the EVD. The rest of the paper is organized as follows. In Section 2, we develop the model. The basic properties (e.g. positivity, boundedness and global existence of the solution) are given in Section 3. Sections 4 and 5 deal with the existence and stability of equilibrium points for the sub-model and the full model, respectively. A dynamically consistent NSFD scheme is presented in Section 6, which also includes numerical simulations that support the theory. Concluding remarks on how our findings fit in the literature and on potential extensions are given in Section 7.

Model formulation
To address the fundamental research question of explaining the recurrence and persistence of EVD outbreaks in Africa from the consumption of contaminated bush meat, the funeral practices and the environmental contamination, we make the following assumptions: (1) Deceased human individuals can still infect during funerals. This assumption is motivated by the African practices (e.g. washing of deceased individuals) during burial ceremonies. (2) Infectious and deceased human individuals shed the environment, especially through their urine and stool. This happens in regions with poor sanitary facilities and/or in regions where people do not observe appropriate hygienic practices. (3) Human individuals can be infected through contaminated environment. This central assumption is supported by several works. Firstly, we have the work [31], which established the survival of filoviruses in liquids, surfaces and glasses. Secondly the findings in [7] clearly demonstrated the persistence of Ebola virus in the environment while the environmental contamination was evidenced in [38]. More importantly, it was reported in [9] that human epidemics took off not only by direct contact via bodily fluids but also by indirect contact with contaminated surfaces. In particular, it was already observed in [17] that an individual contracted the Ebola virus in Uganda after using a blanket previously belonging to a positive case. The above-quoted references show further that epidemiological data support the possibility of nosocomial transmission through indirect contact with contaminated surfaces and fomites. (4) We assume that the consumption of (contaminated) bush meat may result in the provision of Ebola viruses in the environment at a constant rate. This assumption makes sense since, in Africa, and particularly in the regions that were affected by Ebola outbreaks, people live close to the rain-forests, hunt bats and monkeys for meat and go to the forest to harvest fruits for food [23,24], (5) There is a permanent disease-induced immunity. This assumption is motivated by the fact that it has never been reported that an individual caught the EVD for the second time. (6) We assume that there is a vital dynamics. Indeed, some of the Ebola outbreaks have lasted more two years (for instance the Western Africa outbreak). Thus, during this relatively long period of time, there might be new births or inflow of susceptible individuals from other/surrounded places as well as natural deaths, which allow a demographic process to take place, as studied in [1,19].
Based on the above-mentioned assumptions and motivated by the works in [7,17,23,24,31,38], we develop a new deterministic model as follows: The susceptible human population is replenished by a constant recruitment at rate π . Susceptible individuals S may acquire infection after effective contacts β 1 with infectious and β 2 with deceased human individuals. They can also catch the infection through contact with a contaminated environment at rate λ. Infectious individuals I experience an additional death due to the disease at rate δ and they are recovered at rate γ . Deceased human individuals can be buried directly during funerals at rate b. Susceptible, infectious and recovered individuals die naturally at rate μ.
As far as the assumption ( 4) above is concerned, the environment is contaminated by the Ebola virus at a constant rate σ . Actually the inflow σ can be regarded as a kind of 'black box', which encompasses the provision of Ebola viruses in the environment by all means. In this work, we focus on such a provision via the consumption of contaminated bush meat and fruit bats. Although the recruitment rate σ can be modelled explicitly by incorporating the dynamics of fruit bats and wildlife, we postpone this aspect to another project with a more general model. Moreover, infectious and deceased human individuals can shed the environment at rates ξ and α, respectively. The model parameters are summarized in Table 2.
Following the flow chart in Figure 1, the dynamics of the model is governed by the system of ordinary differential equations given by Figure 1. Ebola transmission diagram: black arrows represent the classical transfers in, out or inbetween for different compartments; the lines without arrows represent the interactions leading to infection. The arrow from S to I represents both the infection and the in-between transfer. The lighter arrows from I to P and from D to P represent the shedding.
In principle, the mass action law is used for mathematical simplicity. However, this mass action principle has some biological relevance for the case under consideration given the correlation between the huge gatherings of people at funerals and the fast spread of EVD in the region. We stress that the proposed model (1) is new and extends the works [1,2,10,15,16,19,21,22,30,32,35,37], as discussed in the introduction and also at the beginning of this section. In particular, in considering infection through contact with deceased humans (during funerals), our model differs from [30], where an epidemiological class of death individuals was considered, with nodisease-induced death rate, in that, we do not consider deceased individuals as an explicit compartment. However, these individuals are assumed to remain infectious until they are buried.
System (1) is obviously appended with initial conditions Adding the first, the second and the third equations of (1), we obtain the conservation law where H = S+I+R is the total active human population, i.e. individuals that are alive. Throughout this work, we make the natural additional assumption that the burial rate b is less than or equal to the overall death rate μ + δ. If this condition is not met, then the deceased human individuals completely disappear from the community and the consideration of the compartment (D) in the model is irrelevant.

Well-posedness of the model and equilibria
In this section, we prove that model (1) is well-posed. This is done in three steps below.

Well-posedness
Proposition 3.1: Assume that model (1) has a global solution corresponding to non-negative initial conditions. Then the solution is non-negative at all time.
The first equation of model (1) can be written as

This is a linear first-order equation in S which has solution
Hence S(t) ≥ 0, ∀t ≥ 0. Regarding the non-negativity of the remaining variables, we consider the subsystem which can be rewritten in the matrix form where We note that M is a Metzler matrix (i.e. with non-negative off-diagonal entries) in view of the already established non-negativity of S. Thus, Equation (4) is a monotone system [34]. Therefore, R 4 + is invariant under the flow of system (4). This completes the proof of the proposition.
We can now state and prove the following proposition, which guarantees the boundedness of the solutions of system (1).

Proposition 3.2:
Assume that the initial conditions for system (1) satisfy Then, whenever the solution exists on an interval J, it satisfies the following a priori bounds: Proof: Since, I(t) ≥ 0, we have from Equation (2) that Application of the Gronwall inequality yields from which Consequently, I(t) ≤ H m . Replacing this in the fourth equation of (1) gives from where another application of the Gronwall inequality leads to The boundedness of P(t) is proved similarly.
Combining Propositions 3.1 and 3.2 together 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).

Equilibria
In this subsection, we investigate the existence of equilibrium points of model (1). The constant inflow σ of the Ebola viruses from the environment plays a significant role in the study. This is at first seen from the obvious fact that there is no disease-free equilibrium if σ is positive.
π − (β 1 I * + β 2 D * + λP * )S * − μS * = 0, Adding the first and the second equations in (5) and using the other equations, we have Putting the expressions in Equation (6) into the second equation of (5), we obtain after lengthy algebraic calculation, the quadratic equation in I * where In what follows, it is convenient to write A 1 (σ ), the only coefficient that can have different signs, in the form where Hereafter, thanks to Equation (7) and depending on the value of the virus recruitment rate σ , the existence, the number and the bounds of the equilibrium points are discussed.

Model with Ebola virus-free environment
Throughout this section, we assume that there is no recruitment/provision of Ebola virus: σ = 0.

Existence of equilibria
Since A 0 (0) = 0 in this section, I * = 0 (component of the disease-free equilibrium E 0 ) is always a root of the quadratic Equation (7), while a unique positive root of Equation (7) exists if and only if A 1 (0) > 0. Notice from Equation (9) that A 1 (0) > 0 if and only if R 0 > 1. Consequently, using Equations (6) and (10), the unique endemic equilibrium for Equation (1) is given by The following theorem gives the number (depending on the range of R 0 ) of equilibria for model (1).

Theorem 4.1:
(a) The model (1) always has a disease-free equilibrium given by Equation (11).
Notice that the expression of R 0 above can be recovered using the method in [36]. Moreover, following the interpretation in [8], R 0 is the sum of three infection contributions: , the contribution of infectious humans I; , the contribution resulting from manipulation of infected corpses D; • R P 0 = πλ(bξ + αμ + αδ)/bημ(μ + δ + γ ), the contribution due to the environmental contamination by the virus P.
Remark 1: Assume that we use the standard incidence framework or keep the mass action framework with the initial susceptible population set to be unity. If the environmental class (P) and the vital dynamics are neglected, then the model (1) reduces toan SIRD-type model [15,21,30]. The basic reproduction number of such a model is Using the same parameter values as in [1,2,19,21,30], we observe that this R SIRD 0 underestimates the basic reproduction numbers obtained in the literature as displayed in Table 3. Thus, the importance of the environmental contamination on the dynamics of EVD is further reinforced in order to obtain more accurate basic reproduction numbers, an aspect that is however not the main focus of this work.
Note also that in the absence of the modelling assumption (6) (i.e. there no demographic process), our model falls into the group of 'no vital dynamics models'. We then automatically loose the so-called environmental virus-free endemic equilibrium I * . Therefore, the theoretical comparison of the disease at the endemic level (which is another main aim of this work) is not possible.

Stability analysis of equilibria
System (1) is biologically meaningful in the compact and invariant region defined by The local stability of the disease-free equilibrium when R 0 ≤ 1 and its instability whenever R 0 > 1 follow from the well-known result in [5,36]. For the global stability, we have the following result. Using the fact that π = μS 0 , the Lie derivative of L 0 in the direction of the vector field given by the right-hand side of Equation (1) iṡ Choose c 1 and c 2 such that Thus, With this in mind,L 0 becomeṡ Since it is easy to see that the largest invariant subset contained in the set is the disease-free equilibrium E 0 , we conclude by LaSalle's Invariance Principle [20].

Proposition 4.3: The endemic equilibrium E
The proof of this result is postponed to Section 5, where the proof of the LAS of the unique equilibrium E # (σ > 0) is given. Proposition 4.3 will follow by letting σ = 0.

Theorem 4.4:
In the absence of shedding (α = 0) or manipulation of deceased human individuals before burial (ξ = 0) the endemic equilibrium E * exists and is GAS whenever R 0 > 1.

Proof:
We consider a Volterra-type Lyapunov function where a 1 , a 2 , a 3 and a 4 are four non-negative constants to be determined shortly.
The derivativeL 1 of L 1 along the trajectories of model (1) iṡ Since E * is an equilibrium point, the following relations hold: (μ + δ + γ )I * = β 1 S * I * + β 2 S * D * + λS * P * , Using the relations in Equation (12),L 1 becomeṡ After the expansion and grouping of the expression above, we havė Choose a 1 , a 2 , a 3 and a 4 such that the expressions in the brackets vanish. That is Fixing a 1 = bη and exploiting Equation (12) to solve Equation (14) for a 2 , a 3 , a 4 yields Replacing Equation (15) in (13) giveṡ Grouping some terms in the expression above yieldṡ Using Equations (12) and (15), one can see that Putting this in Equation (16) and performing some algebraic manipulations, we obtaiṅ Henceforth, Case 1 if α = 0, then we havė With this in mind, Equation (17) can be rewritten aṡ In both cases, using the arithmetic-geometric mean inequality, one conclude thatL 1 ≤ 0. Thus, L 1 is indeed a Lyapunov function. Furthermore, is the endemic equilibrium E * . We conclude by LaSalle's Invariance Principle that E * is GAS [20].

The full model with σ > 0
The terminology 'full model' refers to the fact that there is a recruitment/provision σ of Ebola viruses. Thus, we assume throughout this section that σ > 0.

Existence and stability of equilibria
The number and local asymptotic stability of equilibria are investigated in the following result: Theorem 5.1: The system (1) has a unique endemic equilibrium point denoted by E # = (S # , I # , R # , D # , P # ), where I # is the unique positive root of the quadratic Equation (7). The equilibrium point E # is LAS.
Proof: Since the coefficients A 2 , A 1 and A 0 are positive, it follows from Descarte's sign rule that Equation (7) has a unique positive root From this relation and Equation (6), we obtain the other coordinates of E # , namely: The local asymptotic stability of the unique steady-state E # is established using the Routh-Hurwitz criterion, as detailed in the appendix.
The stability in Theorem 5.1 can be partially improved as follows:

Proposition 5.2:
In the absence of shedding and manipulation of deceased human individuals before burial, i.e. ξ = α = 0, the unique endemic equilibrium E # is GAS.

Proof:
We consider the Volterra-type Lyapunov candidate function Q = Q(S, I, D, R, P) where the constants k 1 and k 2 will be determined shortly. Since E # is an equilibrium, the following relations hold: Using Equation (19) in the expression ofQ, the derivative of Q along the trajectories, it can be shown as in the proof of Theorem 4.4 thaṫ Now, we choose k 1 and k 2 such that or equivalently With the relations (19) and (20) in mind, Using the arithmetic-geometric mean inequality, we haveQ ≤ 0, which means that Q is indeed a Lyapunov function. Furthermore, The conclusion follows from LaSalle's Invariance Principle, which gives the global asymptotic stability of E # [20].

Remark 2:
We note that the model developed in this paper does not take into account any Ebola control strategy. Furthermore, the epidemic in the Western Africa has persisted for more than two years. Also, the estimated basic reproduction numbers are generally greater than the unity, as reported in the introduction. In view of these three facts, it is reasonable to expect that the epidemic, as modelled here, establishes itself at an endemic level, as shown in Theorem 4.4 and Proposition 5.2. These findings should not be regarded as inconsistent with the reality observed on the field (where the disease outbreak in some places fade out after some period of times). Typically, the relative quick clearance of the EVD is actually due to some adequate control measures put in place by healthcare authorities, an aspect which, we repeat, is not incorporated in the current model.
The epidemiological implication of the consumption of contaminated bush meat (i.e. σ > 0) is an increase in the endemic level of the disease, as described in the next result.

Proof: We show that
Then, on the one hand, On the other hand, one can easily show that Thus, if But, direct calculations show that This proves the theorem.

Remark 3:
As mentioned earlier, Theorem 5.3 provides interesting answers to the main research question considered in this paper. To be more specific, we assume that human individuals consume contaminated bush meat and fruit bats while there is an outbreak of EVD. Then the endemic level of the disease increases with the said consumption to the extent that the number of infected individuals is larger than the corresponding number of infected individuals when there is no such consumption: I # > I * . Therefore, a natural measure of lowering down the severity of the disease is to educate people not to eat bush meat and/or fruit bats. Looking at this in conjunction with Proposition 5.2 and in accordance with the comment at the end of Remark 2, the disease can go to extinction if the educational programme is broadened against funeral practices and environmental contamination in order to have ξ = α = 0. This educational control programme paid off during the 1976 DRC (Zaïre) Ebola outbreak: the disease was cleared out of the afflicted population.

NSFD method
In this section, we design an NSFD scheme [3] that replicates the dynamics of the continuous model (1).

Construction of the scheme
Let X n = (S n , I n , D n , R n , P n ) T denote an approximation of X(t n ) where t n = n t, with n ∈ N, h = t > 0 being the step size. We propose the NSFD scheme for the model (1), where The discrete method (21) is indeed an NSFD scheme because it is constructed according to Mickens' rules [26] formalized as follows in [3]: Rule 1. The standard denominator h = t of the discrete derivatives is replaced by the complex denominator function in Equation (22) which satisfies the asymptotic relation Note that the denominator function φ is expected to better capture the dynamics of the continuous model through the presence of the underlying parameters μ, δ and γ . In fact, exact schemes for a wide range of dynamical systems involve such complex denominator functions [25,33].
Rule 2. Nonlinear terms in the right-hand side of Equation (1) are approximated in a non-local way. For instance, we have P(t n )S(t n ) P n S n+1 insteadof P(t n )S(t n ) P n S n .

Analysis of the scheme
Theorem 6.1: The NSFD scheme (21) is a dynamical system on the biological feasible domain K of the continuous model (1).
Proof: First, we prove the positivity of the scheme (21). It is easy to show that the NFSD scheme (21) takes the explicit form , where we set Thus S n+1 ≥ 0, I n+1 ≥ 0, R n+1 ≥ 0, D n+1 ≥ 0 and P n+1 ≥ 0, whenever the discrete variables are non-negative at the previous iteration. It remains to prove the positive invariance of K. Adding the first, the second and the third equations in (21), one has Therefore, The a priori bounds for D n+1 and P n+1 follow readily from the fact that D n+1 and I n+1 are less than or equal to H n+1 . This completes the proof.
Following the approach in [3] (Theorem 15) and in [11] (Theorem 2), it can be shown that the NSFD scheme preserves the properties of the continuous model (1) given in Proposition 4.3 and in Theorem 5.1 in the precise manner described below.

Theorem 6.2:
The discrete scheme (21) preserves the equilibrium points of the continuous model (1). That is, on the one hand, the NSFD scheme (21) preserves the equilibrium points of the continuous model (1) in the sense that the only fixed points of the scheme (21) are either the endemic equilibrium point of the continuous model (1) when σ > 0 or its disease-free and endemic equilibria when σ = 0. On the other hand, the fixed points have the same stability properties as the equilibrium points.

Theorem 6.3: (1)
The disease-free fixed point (resp. the endemic fixed point ) of the NSFD scheme (21) for the model without recruitment/provision of Ebola viruses is GAS whenever R 0 ≤ 1 (resp. whenever R 0 > 1). (2) The endemic fixed-point of the NSFD scheme (21) for the full model is GAS.
Proof: Let (X n ) ∈ R 5 + be the bounded sequence defined by the NFSD scheme (21). We want to prove that X n tends to X * , where X * is any of the fixed point in Theorem 6.3. By Bolzano Weierstrass theorem, there exists a subsequence (X n k ) of (X n ) that converge to some Y * as k → +∞.
By the assumption made above and the structure of the NSFD scheme (21), Y * = X * is necessarily either the unique disease-free fixed-point E 0 (whenever R 0 ≤ 1) or the unique endemic fixed-point E * or the unique endemic E # , which is LAS thanks to Theorem 6.2. Therefore, there exists θ > 0 such that for an initial condition X 0 satisfying Let X 0 be an arbitrary initial condition. As lim k→+∞ X n k = X * , there exists an integer k 0 such that X n k 0 − X * ≤ θ .
This proves that X * is GAS.
We conclude this section by providing some numerical simulations. The parameters displayed in Table 4 are mostly taken from recent works on the Western Africa Ebola outbreaks [2,15,35]. The simulations are performed using the NSFD scheme (21) and coded with MatLab. In contrast to Figure 2(a) where Ode45 exhibits negative solutions, Figure 2(b) illustrates the power of the NSFD scheme (21) to produce positive solutions   for any value of the step size h (h = 4). Figure 3(a) shows the global asymptotic stability of the disease-free equilibrium E 0 as established in Theorem 4.2. Figure 3(b) supports the global asymptotic stability of the endemic equilibrium E * in Theorem 4.4.
In accordance with Theorems 5.2 and 6.3, the global asymptotic stability of the equilibrium E # is displayed in Figure 4(a). Figure 4(b) compares I * and I # as demonstrated in Theorem 5.3. For Figure 4(b), we choose the same set of parameters as in Figure 3(b), except that σ = 0 (blue) and σ = 0.6 > 0 (red).

Conclusion
The transmission dynamics of the Central and Western Africa Ebola outbreaks has been relatively more studied in the direct (human-to-human) route [1,2,10,15,16,19,21,22,30,32,35,37] than in the indirect environmental (environment-to-human-to-environment) route. In this paper, we have focused on the latter route, which is one of the main characteristics of the EVD, specifically in Western Africa [7,9,31,32,38]. Rather than dealing only with the direct transmission, we formulated the main research question in the following comprehensive manner in order to incorporate the unavoidable indirect transmission route: Can the consumption of contaminated bush meat, the funeral practices, and the environmental contamination explain the recurrence and persistence ofEVD outbreaks in Africa? More precisely, we enriched the classical SIR model with two additional compartments of deceased individuals who undergo funeral practices and free-living viruses with vital dynamics. In this double setting of direct and indirect transmissions, we made major contributions in two directions.
From the theoretical point of view, we showed in the following precise manner that the severity of the disease increases with the recruitment/provision of Ebola viruses and the disease dies out in the absence of such recruitment/provision as well as in the absence of shedding and manipulation of deceased individuals: (i) The full model has only one endemic equilibrium, which is locally asymptotically stable (LAS) while in the absence of shedding or manipulation of deceased individuals, this equilibrium is GAS. (ii) For the model without recruitment/provision of Ebola viruses, the disease-free equilibrium is GAS when the threshold quantity is less than or equal to unity. When the said threshold is greater than one, the endemic equilibrium is LAS in the generic case and GAS in the absence of shedding and manipulation of deceased individuals.
(iii) At the endemic level, for both cases in items (i) and (ii) above, the number of infected individuals when there is provision of Ebola viruses is larger than the corresponding number of infected individuals in the absence of such provision. (iv) Along the lines of the above-mentioned findings, there is a natural control measure to low down the severity of the disease and even to eradicate it. That is to educate people not to eat bush meat and/or fruit bats, to avoid unsafe funeral or burial ritual practices. This educational programme paid off during the 1976 DRC (Zaïre) Ebola outbreak where the disease was cleared out of the afflicted population.
From the numerical analysis and computational point of view, we introduced the NSFD approach which to the authors' best knowledge has never been applied to the modelling of EVD. We proved analytically and illustrated by numerical simulations that our NSFD scheme is dynamically consistent with respect to the following properties of the continuous model: positivity and boundedness of solutions, local and global stability of equilibria.
As the modelling of EVD is not sufficiently developed, this work offers many opportunities for improvements and extensions. Currently, we are focusing on the following aspects: (a) the use of other incidence functions (standard for the direct transmission, and/or Holling-type functions for the indirect transmission); (b) the extension of the model to account for the complex ecology of the EVD transmission; (c) the multi-species setting to account for the extreme case where a region is attacked by more than one Ebola virus species; (d) the effective modelling of the source of the viruses in the environment by incorporating the dynamics of wildlife and fruit bats; (e) the incorporation of the human behaviour in the model via educational campaigns and media broadcasting; (f) the modelling of some optimal control strategies, such as early detection, isolation and quarantine, treatment, and sterilization techniques; and (h) the incorporation of patches to account for the circulation of the disease in many counties and/or countries as it is the case in Western Africa [6].

Acknowledgments
This work originates from group projects assigned to participants in the 2nd Joint-UNISA-UP