Stability analysis of a compartmental SEIHRD model for the Ebola virus disease

In this work, we perform a stability analysis of a compartmental SEIHRD model. This model is a simplified version of a previous approach. In this previous work, we proposed an original deterministic spatial-temporal model, called BeCoDiS (Between-Countries Disease Spread), to study the evolution of human diseases within and between countries. This model was validated by considering data from the 2014-16 West African Ebola Virus Disease epidemic. Here, considering some simplification assumptions in Be-CODIS, our goal is to study the equilibria of the model and their stability using the basic reproduction ratio as a threshold parameter. Finally, we validate the obtained results by considering some numeriDiène Ngom Département de Mathématiques, Université Assane Seck de Ziguinchor, BP 523 Ziguinchor, Senegal, UMI 2019-IRD & UMMUSCO-UGB e-mail: dngom@univ-zig.sn Benjamin Ivorra Instituto de Matemática Interdisciplinar, Departamento de Análisis Matemático y Matemática Aplicada and MOMAT Research Group, Complutense University of Madrid, Plaza de Ciencias, 3, 28040, Madrid, Spain e-mail: ivorra@mat.ucm.es Ángel M. Ramos Instituto de Matemática Interdisciplinar, Departamento de Análisis Matemático y Matemática Aplicada and MOMAT Research Group, Complutense University of Madrid, Plaza de Ciencias, 3, 28040, Madrid, Spain e-mail: angel@mat.ucm.es Citation: Diène Ngom, Benjamin Ivorra, Ángel M. Ramos, Stability analysis of a compartmental SEIHRD model for the Ebola virus disease, in R. Anguelov, M. Lachowicz (Editors), Mathematical Methods and Models in Biosciences, Biomath Forum, Sofia, 2018, pp. 44-56, http://dx.doi.org/10.11145/texts.2017.12.165 Copyright: c © 2018 Ngom et al. This article is distributed under the terms of the Creative Commons Attribution License (CC BY 4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.

cal experiments based on data from the 2014-16 West African Ebola Virus Disease epidemic.

Introduction
Modeling and simulation are important decision tools that can be used to control human and/or animal diseases [1,13,21].In a previous work (see [14]), we have presented a spatial-temporal epidemiological model, called Be-CoDiS (Between-Countries Disease Spread), for the study of the spread of human diseases between and within countries.This model is an adaptation of Be-FAST (Between Farm Animal Spatial Transmission), which simulates the spread of animal diseases between farms [13,21,19,20,18].More precisely, Be-CoDiS is based on the combination of a deterministic Individual-Based model (where countries are considered as individuals) [8], simulating the between-country interactions (here, migratory flux) and disease spread, with a deterministic compartmental model [5] (a system of ordinary differential equations), simulating the within-country disease spread.This model also considers dynamic parameters to tackle the application of possible sanitary control measures.Be-CoDiS was validated by considering the case of the 2014-16 West African Ebola Virus Disease (EVD) epidemic [11,6].EVD is a human and primates virus disease that causes a high mortality rate in affected population (between 50% and 90%) [10].
In Ref. [14], we have made a detailed description of the model and have presented some numerical results.In this previous article, no theoretical analysis of the mathematical properties of Be-CoDiS was performed.Thus, here, our goal is to propose a novel analysis regarding the mathematical stability of a simplified version of Be-CoDis, focusing only on one country with constant parameters.This simplified approach is based on a compartmental SEIHRD model.To this aim, we analyze the behavior of the equilibrium states of this model.In particular, we estimate the disease basic reproduction ratio , denoted by R 0 , according to the model parameters (see [3] for more details).From a mathematical point of view, it is generally observed that if R 0 > 1 then the epidemic becomes endemic, whereas if R 0 ≤ 1 then the epidemic disappears [1].Finally, we validate the obtained theoretical results by considering numerical experiments based on data from the 2014-16 West African Ebola virus epidemic.We note that similar works have been proposed in the literature for other Ebola outbreaks and models (see, e.g., [7,17]).However, as the structure of those models was different from the one used here and the evolution of the considered epidemics was not similar to the 2014-16 one (which is considered as exceptional, due to its magnitude), the study proposed here is of high interest to understand the possible changes in the dynamic of Ebola epidemics.
This work is organized as follows.In Section 2, we formulate the considered SEIHRD model.In Section 3, we study its equilibrium states.In Section 4, considering data from the 2014-16 West African Ebola virus epidemic, we validate and illustrate the theoretical results with numerical experiments.

Mathematical modelling
We consider a disease with the following states for people (see [14]): • Susceptible (denoted by S): The person is not infected by the disease pathogen.
• Infected (denoted by E): The person is infected by the disease pathogen but cannot infect other people and has no visible clinical signs (e.g., fever, hemorrhages, etc.).After an incubation period, the person passes to the Infectious state.For the sake of simplicity, we assume that the population size in the considered country is constant and equal to N ∈ N (i.e., death flows are compensated by birth flows entering the susceptible state).This hypothesis is reasonable as, due to the size of the population in a country and the time scale of the study (generally lower than . five years) considered here [12].A validation of this assumption can be found in [14].Furthermore, to simplify the notations, we consider that S, E, I, H, R and D denote the ratio of people in each state in the considered country (rather than the total number of people).A diagram of this model for one country is shown in Figure 1.
Under these assumptions, the evolution of the epidemic, is modeled by where

Stability Analysis
Let It is straightforward to see that the set Ω is positively invariant for System (1) (see [4]).
For the study of the stability properties of System (1), we will compute and use the basic reproduction ratio.
From a mathematical point of view, the value of R 0 associated to the epidemiological compartmental model (1) can be computed as the spectral radius of the so-called next generation matrix (see [23] for more details and notations).
We introduce the following matrix formulation of System (1).Let P = (X,Y ) T , with X = (E, I, H, D) T and Y = (S, R) T .System (1) can be rewritten as where and g(S) defined by Then, we state the following theorem.
Theorem 1.The basic reproduction ratio of System (1) (or, equivalently (2)) is given by Proof.To compute the basic reproduction ratio of the considered system, we apply the Next Generation Matrix methodology briefly described above (see [23]).
To do so, we first determine the disease free equilibrium points of System (1) in Ω , i.e., the points of the form P = (X,Y ) T , with X = (0, 0, 0, 0) T and Y = (a, b) T such that a ≥ 0, b ≥ 0, a+b = 1, F (S) − Ṽ X = 0 and g(S)P = 0.It is straightforward to see that P f = (X f ,Y f ) T , with X f = (0, 0, 0, 0) T and Y f = (1, 0) T , is the unique admissible disease free equilibrium point of this system in Ω (we note that all the hypothesis required for using this methodology, detailed in [23], are satisfied).
We consider , and obtain that F = F (1) and V = Ṽ .After computation, we observe that F and V −1 are non-negative matrices.Thus, following the technique used to compute the value of the basic reproduction ratio presented in [9], we obtain after some computations (using the software Maple 16) that the value of R 0 is given by R 0 = ρ(FV −1 ) which leads to Using this basic reproduction ratio R 0 , we have the following stability results.
After some computation, using Maple 16, we obtain that P f and P e , defined previously, are the equilibrium points of this system in Ω .
Let us assume that R 0 ≤ 1: We use a method developed in Ref. [22] to determine a Lyapunov function for the disease free equilibrium P f .Considering this goal, the first line of System (2) is rewritten as where F and V are defined previously and Symbolic computations done with Maple 16 show that R 0 is an eigenvalue of matrix V −1 F. Furthermore, is a left eigenvector, with eigenvalue R 0 , of the matrix V −1 F. Let L f : R 4 × R 2 → R be given by By calculation, we obtain that L f is non negative in the set [0, 1] 6 , L f (P f ) = 0 and is non negative for all (X,Y ) ∈ Ω .Additionally, since R 0 ≤ 1 and the coordinates of wV −1 and f (X,Y ) are non negative for all (X,Y ) T ∈ Ω , then Lf (X,Y ) ≤ 0 for all (X,Y ) ∈ Ω .Hence, L f is a Lyapunov function of System (2) at the equilibrium state P f and, thus, the equilibrium P f is globally stable since the hypotheses of the Lyapunov Theorem are satisfied in Ω (Lyapunov 1892; see, for instance, page 5 in Ref. [15]) .Moreover, to show that P f is globally asymptotically stable, we use the LaSalle's Invariance Principle (see, for instance, Theorem 2 in Ref. [16]).Let We obtain that Lf (X,Y ) = 0 if and only if (R 0 − 1)wX = 0 and wV −1 f (X,Y ) = 0.This implies that S = 1 or I = H = D = 0. Thus, Γ f = {(1, 0, 0, 0, 0, 0 Therefore, dS(t) dt As δ > 0, this leads to Let Γ f ,0 be the largest invariant set of System (1) in Γ f .In {(S, E, I, H, R, D) T ∈ Ω /E = I = H = D = 0}, System (3) can be rewritten as and the solution of this system is given by S(t) = S(0) • e −µt + 1 > 1, which is absurd as 0 < S(t) ≤ 1. Hence Γ f ,0 is reduce to P f = (1, 0, 0, 0, 0, 0) T .Due to the Lasalle's principle, we conclude that P f is globally and asymptotically stable.
Focusing on the equilibrium point P e , if R 0 ≤ 1 then P e is not admissible in Ω .
Let us assume that R 0 > 1: System (2) satisfies the hypothesis of Theorem 2 in Ref. [23] and, thus, the equilibrium point P f becomes unstable when R 0 > 1.
We now focus on the study of the equilibrium state P e , which corresponds to the endemic equilibrium.Since S(t) + E(t) + I(t) + H(t) + R(t) + D(t) = 1, we can remove the second equation of System (1).In this case, the linearized version of this system at point P r e = (S e , I e , H e , R e , D e ) T can be written as where Z = ( Ŝ, Ŵ ) T , with Ŝ = S−S e , Ŵ = W −(I e , H e , R e , D e ) T and W = (I, H, R, D) T , and M (P r e ) a matrix associated to System (1).

Numerical Experiments
In this section, in order to validate and illustrate the interest of the theoretical results obtained previously, we present some numerical experiments based on data from the 2014-2016 West African EVD epidemics [11,6].
To study some representative numerical examples, we consider two set of parameters, denoted by Set 1 and Set 2, detailed in Table 1.Set 1 and Set 2 are within the range of values proposed in [14] for the 2014-2016 West African EVD case and correspond to cases associated with basic reproduction ratios of 0.3291 and 1.3910, respectively.
The initial conditions are set to S(0) = 0.999, E(0) = 0.001 and all other ratios set to 0. The model is discretized by considering an explicit Euler scheme with a step size of 0.1 day.The simulation is stopped after a maximum number of 3650 days; or if the evolution of people in state S from one iteration to other is lower than Taking into account those parameters and numerical methods, we perform the following two experiments: • Experiment 1: We consider the set of parameters Set 1.The evolution of the ratio of contaminated people is presented in Figure 2. In this case, this ratio is decreasing.The simulation stops after 91 days due to the low ratio of contaminated people.
• Experiment 2: We consider the set of parameters Set 2. The evolution of the ratios of contaminated and safe people (i.e., people either in the state S or R) are shown in Figure 3.We can observe that the epidemic reaches an endemic equilibrium of 13% of contaminated people in the population.The simulation stops after 1149 days due to the stabilization of the numerical solution.The endemic equilibrium satisfy S e = 0.76534, E e = 0.04725, I e = 0.010609, H e = 0.079172, R e = 0.097633 and D e = 0, which is numerically close from the theoretical results presented previously.

Conclusion
In this paper, we have performed an analysis of the equilibrium states of a compartmental SEIHRD model, corresponding to a simplified versions of the Be-CoDiS model proposed in [14].First, we have estimated the basic reproduction ratio (denoted by R 0 ) of this model.In particular, we have obtained in Theorem 2 an analytical expression of R 0 according to the model parameters.Additionally, we have proven that if R 0 ≤ 1, then the disease free equilibrium is globally and asymptotically stable, which is a desirable biological situation because the epidemic will disappear.When R 0 > 1, we shown that the disease free equilibrium is unstable and the endemic equilibrium is locally stable.This leads to the persistence of the epidemic in the considered population.Those results have been validated with representative numerical experiments.One of the main interest of those results is to propose a decision tool in the case of future outbreaks.Indeed, as said previously, the behavior of the evolution of a possible Ebola epidemic depends on the value of R 0 .This basic reproduction ratio was expressed in function of the model parameters.By performing a sensitivity analysis of those parameters on R 0 , we can determine which ones have the highest influence on the dynamic of the epidemic and allocate the sanitary resources to control, if possible, those parameters.
In future works, as we aims to study the spread of human diseases between countries, we will complete the analysis of this simplified SEIHRD model by considering several countries and a migration flow between them.Additionally, we will consider whether the dimensionality of the model could be reduced (See, e.g., [2]).

Fig. 1
Fig. 1 Diagram of the simplified model for one country considered in Section 2

Fig. 2
Fig. 2 Evolution of the ratio of contaminated people simulated during Experiment 1 presented in Section 4.

Fig. 3
Fig. 3 Evolution of the ratios of contaminated and safe people simulated during Experiment 2 presented in Section 4.
The person is hospitalized and can still infect other people.At the end of this state, the person can pass either to the Recovered state or to the Dead state.We point out that state H does not include hospitalized people which cannot infect other people any more.This last category of people is included in the Recovered state explained below.•Dead (denoted by D): The person has not survived the disease.The cadavers of infected people can infect other people until they are buried.After a fixed average time, the body is buried and can no longer infect other people.• Recovered (denoted by R): The person has survived the disease, is no longer infectious and develop a natural immunity to the disease pathogen.
α, λ and θ denote the transition rates (day −1 ) from a person in state E to state I, from state I to state H, from state H to state R, from state H to state D and from state D to state S, respectively.

Table 1
Values of the parameters in Set 1 and Set 2 used in during the experiments presented in Section 4. The basic reproduction ratio (R 0 ) generated by those values is also reported below.