Conditions for endemicity: qualitative population dynamics in a long-running outbreak of Ebola virus disease

Abstract Ebola virus disease (EVD) struck West Africa in 2013–2016 in an epidemic of unprecedented scope, with over 28000 cases and 11000 fatalities in the affected region. The protracted duration of the outbreak – more than two-and-one-half years of active transmission – raises questions about the persistence of EVD. In this brief paper, we qualitatively examine conditions supporting long-running EVD epidemics via a susceptible – exposed – infectious – recovered – deceased-infectious differential equations model that incorporates births and non disease-related deaths. We define an ‘effective epidemiological population’ to include contagious individuals recently deceased from the disease. Under a constant effective epidemiological population condition, we consider the basic reproductive number and use Lyapunov function arguments to establish conditions in the parameter space supporting an exchange of stability from the disease-free equilibrium to an endemic equilibrium.


Introduction
Ebola virus disease (EVD) is a contagious disease in humans caused by Ebolavirus, a genus in the family Filoviridae. In the period between its first identification in 1976 and 2013, EVD appeared only in isolated outbreaks, accounting for fewer than 2400 cases and 1600 fatalities in total (Centers for Disease Control and Prevention, 2017). In December 2013, EVD struck again in Guinea, the beginning of an epidemic of unprecedented scale and duration (World Health Organization, 2016a). Over the course of two-and-one-half years, three countries in the disease epicenter -Guinea, Liberia and Sierra Leone -witnessed an EVD caseload of over 28000 infections leading to more than 11000 deaths. The West African epidemic earned two somber distinctions: it is the largest and most lethal EVD event to date, as well as the longest running.
The protracted nature of the West African epidemic presents challenges to its analysis. In the study of outbreaks of short duration, mathematical models are typically constructed with the assumption that the contribution of background population dynamics is negligible CONTACT Matthew Glomski matthew.glomski@marist.edu on a brief epidemiological timescale. Omitting these demographic effects, such as births and non disease-related deaths, is an often justifiable and worthwhile simplification. In this paper, we choose instead to incorporate underlying population dynamics in our compartmental model, and qualitatively examine EVD epidemics of long duration.

Background on compartmental models for EVD
Compartmental models describing the transmission of EVD appear in the literature as early as 1996. Astacio et al. (1996), and Chowell, Hengartner, Castillo-Chavez, Fenimore, and Hyman (2004) in their seminal paper, used a susceptible-exposed-infectious-removed (SEIR) model to study EVD outbreaks in the Democratic Republic of the Congo (DRC) and Uganda. The SEIR model is effective in describing disease transmission from living infectious individuals to the susceptible population. Yet there is overwhelming evidence of an additional pathway of EVD infection: one in which transmission also occurs via contact with those recently deceased from the disease (e.g. Hewlett & Amola, 2003;Francesconi et al., 2003;Peters, 2005;Allaranga et al., 2010;Richards et al., 2015). The first model to incorporate a 'deceased-infectious' compartment may be that from Legrand, Grais, Boelle, Valleron, and Flahault (2006). They examined the 1995 DRC and 2000 Uganda EVD outbreaks through an SEIHFR compartmental model, wherein susceptible individuals are at risk via contact with the infected (living) population in the community (I) or hospital (H), or via contact during funerals and/or burials (F) with those recently deceased from the disease. The spirit of this approach has been revisited to describe disease dynamics in the 2013-2016 West African EVD epidemic, and in other filoviral outbreaks. Models from Brozek (2011); Gomes et al. (2014), Rivers, Lofgren, Marathe, Eubank, and Lewis (2014), Camacho et al. (2014), Althaus, Gsteiger, and Low (2014), Lewnard et al. (2014), Pandey et al. (2014) and Weitz and Dushoff (2015) each track the effect of postmortem infection in filoviral epidemics via systems with an explicit deceased-infectious compartment. One excellent review of several of these models is from Chowell and Nishiura (2014).

The case for demographic effects
Previous compartmental models describe an outbreak that is brief with respect to demographic effects, and do not account for the epidemiological impact of births and non disease-related deaths. This is often a justifiable and worthwhile simplification, but one we will not make in our model for the West African EVD epidemic. The demographic timescale is shortened with respect to the epidemiological one by relatively low life expectancies in the affected region. Liberia and Guinea rank 158 and 167, respectively, among 183 countries, whereas Sierra Leone -with a life expectancy at birth of only 50.1 yearsranks last (World Health Organization, 2016b). Simultaneously, these countries feature a high fertility rate, with 4.6 live births per woman in Sierra Leone, 4.7 in Liberia and 5.0 in Guinea. These numbers are high even by comparison to the average for states under the UN designation 'Least Developed Countries', with an average fertility rate of 4.2 (The World Bank, 2014). Considered in the light of a two-and-one-half-year epidemiological event, it is prudent to consider the demographic contributions of births and non disease-related deaths.
The incorporation of background demographic effects in certain other compartmental epidemiological models has been addressed in well-known examples in the literature (e.g. Li and Muldowney, 1995;Hethcote, 2000). Specifically, Lyapunov function arguments are used to qualitatively analyse asymptotic stability in SIS, SIR, SIRS, SEIS and SEIR models in Korobeinikov (2004a), Korobeinikov (2004b) and Korobeinikov & Wake (2002). We adapt some of the techniques used in those papers here, looking to extend these results to a new model.

Constructing an SEIRDi model
We construct a new model for the spread of a long-running epidemic of EVD that incorporates a non-negligible case fatality rate and a route of transmission from the recently deceased to the living. Simultaneously, we include the effect of background demographic dynamics.
We begin with a deceased-infectious (SEIRD-type) model similar to several proposed above (Brozek, 2011;Althaus et al., 2014;Weitz & Dushoff, 2015). Susceptible individuals (S) who become infected move into an exposed class (E). Members of the exposed class are infected, but do not display symptoms of the disease, nor can they transmit it. Exposed individuals move into the infectious class (I) after an exponentially-distributed waiting time with mean κ −1 . Following a second exponentially-distributed waiting time, in this case with mean α −1 , infectious individuals either die from the disease with probability ρ, or move to the recovered class (R) with probability 1 − ρ. We assume recovery confers lifelong immunity. Those infectious individuals who die from the disease move into a deceased-infectious class (Di). Here we use the notation Di, rather than D, to emphasize that this compartment includes only those deceased members of the population that can still transmit infection to the susceptible population. Individuals in the class Di leave the compartment following safe burial or cremation, exponentially distributed with mean waiting time δ −1 . All parameters are assumed to be positive.
The model incorporates no vertical transmission. It has been observed that viral material in the bloodstream of an EVD-infected pregnant woman tends to concentrate in the placenta and amniotic fluid (Mèdecins Sans Frontières, 2015). Live births are rare, and neonatal survival is rarer still (Black, Caluwaerts, & Achar, 2015).
We next consider demographic effects. We assume a constant influx of susceptibles, via births at a rate μ, and subtract from each compartment the same mass-dependent term representing non EVD-related deaths. This approach to modelling background population dynamics is standard in the literature (e.g. Li and Muldowney, 1995;Hethcote, 2000;Korobeinikov, 2004a). See System (1) and (Table 1).  (1).

S(t) proportion of initial effective epidemiological population susceptible E(t)
proportion of initial effective epidemiological population exposed proportion of initial effective epidemiological population (living) infectious R(t) proportion of initial effective epidemiological population recovered with immunity mean time from disease death to safe burial/cremation (days) (1). Susceptible individuals (S) enter the exposed compartment (E) upon infection. Following an incubation period, these individuals next move into the infectious compartment (I). A proportion ρ of the infectious compartment will die from the disease, moving into the deceased-infectious compartment (Di). These deceased individuals remain infectious until safe burial or cremation. Infectious individuals who survive the disease (proportion 1 − ρ) move into the recovered compartment (R), with lifetime immunity. Other flux into and out of each compartment is due to demographic effects.

Figure 1. A schematic diagram depicting EVD transmission as in System
We define N(t) = S(t) + E(t) + I(t) + R(t) + Di(t) to be the 'effective epidemiological population', noting that this count includes deceased-infectious individuals. Depending on assumptions in the parameter space, N(t) may or may not be constant. For simplicity, we choose our five state variables S(t), E(t), I(t), R(t) and Di(t) to represent the fraction of the initial population corresponding to each epidemiological compartment. See also (Figure 1).

Disease-free equilibrium and the basic reproductive number
For the biologically feasible parameter values: μ ≥ 0, κ −1 , α −1 , δ −1 > 0, 0 ≤ ρ ≤ 1, SEIRDi System (1) is well posed. We consider the relationship between background population growth flux rate μ and the inverse mean time from disease death to safe burial δ. The case μ > δ implies dN dt is strictly positive: an a priori assumption of unbounded population growth. We do not believe that this is an accurate reflection of the dynamics behind a long-term epidemic of EVD (or those behind any human epidemic), and therefore, restrict our analysis to the case μ ≤ δ. Under the assumption of biologically feasible parameter values, System (1) admits the disease-free equilibrium Using the next generation method (Diekmann, Heesterbeek, & Roberts, 2009), we calculate the basic reproductive number R 0 . The basic reproductive number is an essential epidemiological threshold quantity, defined as the number of secondary infections caused by the introduction of a single infectious individual into an otherwise wholly susceptible population. Let where T and are the Jacobian matrices of transmissions, and transitions, respectively, about X 0 . With The basic reproductive number R 0 is the dominant eigenvalue of K:

Endemic equilibrium
The system undergoes a bifurcation when the background population flux μ reaches the inverse of the mean deceased-infectious period, i.e. when μ = δ. At precisely this value, the effective epidemiological population N(t) = S(t) + E(t) + I(t) + R(t) + Di(t) = 1 becomes constant. This marks an important threshold, as an endemic equilibrium X * -a persistent non-zero infection state -becomes biologically and mathematically feasible for R 0 greater than unity: (4) The equality of the constants μ and δ are essential to sustain the endemic equilibrium X * . For R 0 greater than unity, the theoretical disease trajectories in (Figure 2) suggest that X * is in equilibrium only if μ = δ. We consider this exchange of stabilities rigorously in the following section.

Stability results
We establish the following theoretical results for the SEIRDi model (1) for a constant effective epidemiological population. Theorem 3.1: is well posed on the region D: (2) the disease-free equilibrium X 0 is globally asymptotically stable on D for R 0 < 1; (3) at R 0 = 1, the stable disease-free equilibrium X 0 coincides with the endemic equilibrium X * ; and (4) the endemic equilibrium X * is globally asymptotically stable on the region D \ X 0 for R 0 > 1.

Proof:
(1) We adapt an argument from Korobeinikov (2004a). Consider the compact set = S, E, I, Di ∈ R 4 : S, E, I, Di ≥ 0, S + E + I + Di ≤ 1 , and the Lyapunov function In R 4 + \ , the derivative is negative, and there is no flux from through S + E + I + Di = 1. The state S(t) = 0 implies dS dt > 0, so solution trajectories with initial conditions in cannot exit via the hyperplane S = 0. For J ∈ {E, I, Di}, we have dJ dt J=0 ≥ 0, so that escape of a solution from through any J = 0 is also impossible; we can conclude region is positively invariant under (1). Since R(t) is given by 1 − S(t) + E(t) + I(t) + Di(t) , an expression for R(t) ∈ [0, 1] can be recovered when the other values are known. It follows that System (1) is well posed on the region D.
(2) Let R 0 ≤ 1. Consider the candidate Lyapunov function The time derivative of W in the context of System (1) is given by At the disease-free equilibrium X 0 , (S) = (1) = 0 and I = 0, so that dW dt is identically zero. In D \ X 0 , the term (S) is strictly negative for S > 0, and the quantity R 0 is negative whenever R 0 = κ βδ + αγρ δ κ + μ α + μ is less than unity.
Thus for R 0 < 1, dW dt is zero at X 0 and negative elsewhere in D for S > 0. For the case S = 0, observe that dS dt reduces to μ > 0, and the Lyapunov argument above can be applied for any t > 0. We conclude that for R 0 < 1, the disease-free equilibrium X 0 is globally asymptotically stable under System (1).
(3) From Equation (4), observe that S * and I * can be rewritten as S * = R −1 0 and Since E * , R * and Di * are zero whenever I * is, it is readily apparent that X * = S * , E * , I * , R * , Di * = 1, 0, 0, 0, 0 = X 0 for R 0 = 1. (4) The global asymptotic stability of the endemic equilibrium for an SEIR model with demographic effects was proven by Li and Muldowney (1995) via an extension of the Poincaré-Bendixson theorem. Later, Korobeinikov & Wake (2002) and Korobeinikov (2004a) proved parallel results using Lyapunov function arguments. We apply some of the techniques of their papers as we complete the proof of Theorem (3.1).
Let R 0 > 1, and defineα = α + μ andκ = κ + μ. We consider the candidate Lyapunov function V : with c = β α + αγρ αδ S * and f = γ δ S * . Melesse and Gumel (2010) point out that similar functions of Goh-Volterra type appear in Lyapunov arguments elsewhere in the literature. The time derivative of V in the context of System (1) is given by At the point X = X * = (S * , E * , I * , R * , Di * ), all derivatives of System (1) vanish, establishing the following four facts: With these substitutions, the derivative reduces to One final simplification gives the following expression for dV dt : The bracketed expressions in lines (5a), (5b), and (5c) can be viewed as differences between the geometric and arithmetic means of sets of non-negative terms. The geometric mean is less than or equal to the arithmetic mean, with equality in D only at S = S * , E = E * , I = I * and Di = Di * . Thus for R 0 > 1, dV dt is zero at X * and negative elsewhere in D \ X * . We conclude that for R 0 > 1, V (t) is a Lyapunov function about the endemic equilibrium for System (1), and this endemic equilibrium is globally asymptotically stable.
Remark 1: Observe that under a zero-fatality condition ρ = 0, System (1) reduces in the limit to the SEIR model with demographic effects (with no vertical transmission) discussed in Korobeinikov (2004a). Remark 2: Under the condition of 100-percent case fatality, (1) can be considered as a type of staged progression model, albeit one in which the second infectious class is comprised entirely of deceased individuals. With the case fatality ratio ρ taken to be unity, the results here are consistent also with the findings in Melesse & Gumel (2010).

Conclusion
The 2013-2016 West African EVD epidemic is notable for its historic scale and long duration. In this paper, we developed a compartmental model to qualitatively describe the dynamics of human-to-human EVD transmission in a long-running epidemic. The proportions of susceptible, exposed, (living) infectious, recovered and deceased-infectious individuals are represented in an SEIRDi system of differential equations. Our model accounts for transmission of infection to the susceptible class from both the living, and deceased, infectious individuals of the effective epidemiological population. Further, it incorporates the background demographic effects of births and non disease-related deaths. Under the assumption of zero population growth in the effective epidemiological population, we showed that the disease-free equilibrium is globally asymptotically stable for R 0 less than unity. For R 0 greater than unity, there exists a unique endemic equilibrium, shown to be globally asymptotically stable on the biologically feasible region away from the disease-free equilibrium.
The model described in this paper may be equally appropriate in describing transmission dynamics of Marburg virus disease (MVD), the only other filoviral disease of humans. Although rarer than Ebola, MVD is similar clinically and epidemiologically (World Health Organization, 2014). Specifically, the high risk of postmortem infection suggests it may be useful to consider our model in the context of MVD outbreaks as well.
The next step in understanding long-running filoviral epidemics is clear. Given the models discussed here and elsewhere in the literature, and given what real-world epidemiological data is known, it should be possible to perform an objective model selection. In the event of another epidemic of such terrible scale and duration, an appropriate model will be essential in crafting the most effective strategies to minimize suffering and loss of life.