Dynamics of Ebola epidemics in West Africa 2014

Introduction

The outbreak of the 2014 Ebola virus epidemic in West Africa, started in late 2013, does not seem to be under control and accurate predictions appear to be extremely difficult. The major reason for this might be due to unstable treatment conditions that provide different reproduction numbers at different periods. However, there are also other challenges related to the mathematical modeling of this epidemic. To address these challenges, several new models have been suggested that show quite different results, we note a few of them published recently^1–8.

In this article we introduce a new model to study the dynamics of the current outbreak by considering the average infectious period as a time-dependent parameter. It is derived from the well studied SIR (Susceptible-Infectious-Recovery) model with time delay (e.g. 9,11), where the decrease in the number of susceptible population in compartment S is the major force stopping epidemics. The susceptible population S is often considered as a whole population. A major drawback of this model, in terms of the current epidemic, is that the population infected constitutes a very small proportion of the total population, a very small decrease in S has almost no effect on compartment I.

We discuss how this drawback could be tackled and introduce a new model that uses only compartment I. This leads to a linear model having some similarities to those models based only on transmission rates from infectious population at different generations (e.g. 6). Our main goal is to fit data by estimating fewer and most influential parameters without considering many other issues like the infectiousness in hospitals and death ceremonies.

This in addition, allows us to have a more robust model with easily interpreted parameters that can be used for more accurate predictions. The main parameter in this model is the average infectious period τ₂ (time from onset to hospitalization) that defines the dynamics of infectious population. This parameter can also be considered as a control parameter in the development of control models dealing with the spread of infection.

We calculate the basic reproduction numbers R₀ for each country (Guinea, Sierra Leone and Liberia) as well as the total Ebola data worldwide. We also provide predictions corresponding to different scenarios by considering different values for τ₂ for future time periods.

Methods

We use the notation I_a(t) for the number of “active” infectious population at time t; it mainly represents the total number of infectious population that are not yet hospitalized. C(t) and D(t) are the cumulative number of infected cases and deaths, respectively. The population density of a country is denoted by 𝒟. This is used in the definition of the infection force of the disease with coefficient β. Moreover, μ stands for the natural death rate of the population, α for the death rate due to disease, and τ₁ for the average latent period (in days) that infected individuals become infectious and τ₂ for the average infectious period (in days).

The main equations of our model are as follows (see Appendix for details):

Here ω is a gamma (cumulative) distribution function (with p.d.f - ω_p) for deaths due to disease⁸; for the values of the parameters see Appendix. We note that there are only three parameters that need to be estimated to fit data for cumulative number of infected and death cases. These parameters are:

Here α and β are continuous variables, τ₂ is a discrete variable with integer values (days).

Basic Reproduction Number - R₀. We calculate the basic reproduction number by considering the stationary states in (1) as follows:

Since the natural death rate μ is close to zero; that is, 1– μ ≈ 1, from (2) we have

Moreover, since $\alpha {\sum }_{i=0}^{{\tau }_2-1}\omega (i)<1$, the reproduction number R₀ ≈ β𝒟τ₂. This means that the reproduction number depends almost linearly on τ₂.

The effective reproduction numbers - R_k, k ≥ 1. The effective reproduction numbers R_k are considered on several consecutive time intervals Δ_k = [t_k, t_k+1], k = 1, 2, …, with corresponding values of τ₂. They are calculated by the same formula as R₀.

Here we make a reasonable assumption that the transmission rate β𝒟 describes the interaction of population (that is, in some sense, related to the local conditions and the life style) and should become relatively stable in the long term for a particular country. Then, the efforts in preventing the spread of infection are mainly observed in the change (decrease) in the value of τ₂.

The infection rate is the key factor defining the dynamics of infectious population. The study⁵ shows that the infection rate is a linearly decreasing function of the total case reported. In 10 the 1995 Ebola outbreak in Congo is considered where the transmission coefficient β is assumed to decrease exponentially due to control interventions. In our model the transmission of infection depends on two parameters β and τ₂. Here β (having slightly a different meaning to what was used in 10) is constant and the dynamics of infectious populations depends on the change in τ₂.

Therefore, to calculate the effective reproduction numbers, we fit data and find the optimal values for α and β, that are constant over the whole period, and optimal values ${\tau }_2^k$ on each interval Δ_k. Then R_k is calculated by formula (2) setting ${\tau }_2^k$.

The sequence of optimal values ${\tau }_2^1$, ${\tau }_2^2$, …, is considered as a method to describe the effectiveness of measures applied for preventing the spread of infection. This sequence very much defines the reproduction numbers on each consecutive time interval and therefore the dynamics of the infected population. It also allows us to consider future scenarios in terms of possible average infectious periods (i.e. times from onset to hospitalization).

Results and discussion

Data were retrieved from the WHO website (http://www.who.int/csr/disease/ebola/situationreports/en/) for the cumulative numbers of clinical cases (confirmed, probable and suspected) collected till 11 November 2014. In all numerical experiments, the second half of the available data for each country is used for fitting the cumulative numbers of infected cases and deaths. The global optimization algorithm DSO in Global And Non-Smooth Optimization (GANSO) library^12,13 is applied for finding optimal values of parameters.

First we consider the whole period of infection in each country and find the best fit in terms of three variables α, β and τ₂ (Problem (DF₁) in Appendix). The results are presented in Table 1. Although from Figure 1 it can be observed that the best fit for Guinea is not as good as for the other cases, these results provide some estimate for the reproduction number R₀ for a whole period of infection till 11-Nov-2014. In all cases (except Guinea), R₀ is around 1.20 and for Guinea - 1.09. We note that the dynamics of infected population is much more complicated (especially in Guinea) which suggests that the reproduction number has been changing since the start of Ebola-2014 in almost all countries. This fact has been studied in 8 in terms of the instantaneous reproduction number over a 4-week sliding windows for each country (see also the next section for different values for τ₂).

Figure 1. The best fits for the cumulative numbers of infected cases and deaths in Guinea, Sierra Leone, Liberia and worldwide by considering three parameters α, β and τ₂ (for the values see Table 1).

The lines represent the best fits, red and black circles represent the data.

The effective reproduction numbers

According to (2), the basic reproduction number is mainly determined by β and τ₂. Since in our model parameter τ₂ takes discrete values (days) it would be interesting to study the change of this parameter over time while keeping β the same for the whole period. This approach makes it possible to consider different scenarios for future developments regarding the change in this parameter and to provide corresponding predictions.

We consider three consequent time intervals Δ_k = [t_k, t_k+1] (k = 1, 2, 3) for each country and find optimal values α, β and ${\tau }_2^k$ (k = 1, 2, 3) (Problem (DF₂) in Appendix). The results are presented in Table 2. The last time point t₄ is 11-Nov-2014. The values of t₁, t₂, t₃ are as follows: 22-March, 23-May and 20-July for Guinea; 27-May, 20-June and 20-August for Sierra Leone; 16-June, 20-July and 07-Sept for Liberia; and 22-March, 23-May and 07-Sept for the total data (World). Each interval Δ_k has its own reproduction number R_k that defines the shape of the best fits presented in Figure 2.

Table 2. Results of best fits: the (effective) reproduction numbers R_k and average infectious period ${\tau }_2^k$ (in days) for different intervals Δ_k, k = 1, 2, 3.

The optimal values for α and β are also provided; they are constant for a whole period.

Country	α	β	R1 (${\tau }_2^1$)	R2 (${\tau }_2^2$)	R3 (${\tau }_2^3$)
Guinea	0.667	0.00527	0.86 (4)	1.25 (6)	1.07 (5)
Sierra L.	0.353	0.00366	1.72 (6)	1.17 (4)	1.17 (4)
Liberia	0.526	0.00688	1.45 (6)	1.23 (5)	0.99 (4)
World	0.489	0.00519	1.04 (4)	1.29 (5)	1.04 (4)

Figure 2. The best fits for the cumulative numbers of infected cases and deaths in Guinea, Sierra Leone, Liberia and worldwide by considering parameters α, β and three consequent time intervals with different values ${\tau }_2^k$, k = 1, 2, 3 (for the values see Table 2).

The lines represent the best fits, red and black circles represent the data.

Figure 3. The cumulative number of infected population according to different scenarios corresponding different values ${\tau }_2^{k+3}$ for time intervals Δ_k.

The starting values of parameters (α, β and ${\tau }_2^k$, k = 1, 2, 3) are in Table 2 (World). The first time interval (Δ₁) is 12/Nov/2014–31/Dec/2014, followed by each next month and the last interval (Δ₅) starts from 1/Apr/2015. The reproduction number for τ₂ = 3 is R = 0.778; it is less than 1 which leads to stabilization. For corresponding reproduction numbers for τ₂ = 4 and 5 see Table 2 (World).

In all cases the effective reproduction number is still greater than 1. In Liberia it shows a decrease from 1.45 to 0.99 and this can be seen in quite a noticeable decrease in the number of cumulative infected cases (Figure 2).

Appendix

Model

The idea behind the model introduced in this paper is related to the SIR model with time delay. Since we are going to implement it on available daily data, a discrete version of this model is considered with the step-size one day. Moreover, since the “recovered” population is not our focus, we will only consider equations related to susceptible (S) and infectious (I) individuals. The most commonly used SIR model¹ in the literature is provided below (see, for example,⁸):

Here λ is the recruitment of the population; μ is the natural death rate of the population; α is the death rate due to disease; γ is the recovery rate; and τ₁ is the latent period that infected individual becomes infectious.

The fraction (1–μ)^τ₁ represents the survival rate of population over the period of [0, τ₁] (in continuous-time case it is equivalent to e^–μτ₁). Below we examine this model in detail and develop an improved model.

Susceptible individuals. Equation (3) describes the dynamics of susceptible individuals S(t). This equation “keeps” the number of infectious individuals I(t) bounded. For example, when the basic reproduction number is greater than 1, there exists⁸ an endemic equilibrium (S*, I*) and S(t) → S*, I(t) → I* as t → ∞. In the case when the birth rate is zero (λ = 0) the relation I(t) → 0 suggests that S(t) → 0.

Thus according to this model the epidemic ends because the number of susceptible individuals S(t) decreases over time and the effective reproduction number (as a function of time) becomes less than 1 at some stage; that is, the number of newly infected population F(S(t), I(t)) decreases thanks to the “enough” decrease in the number of susceptible individuals (while I(t) still increases). This might be applicable to epidemics in early 1900s but it is definitely not applicable to recent ones.

This issue significantly restricts the application of the SIR model for the study of the current Ebola virus epidemic. Below we consider 3 possibilities to overcome this difficulty.

1. The simplest way would be to use a “relatively small” number S(0) for a possible number of susceptible individuals that may become infected. This approach has been implemented in¹ where the total population size in each country (Guinea, Sierra Leone and Liberia) was assumed to be 10⁶ individuals.

2. An interesting (and most reliable in our opinion) approach would be considering “relatively small” number of population S(0) as a variable that needs to be estimated. We have implemented this approach and the results show that currently available curve/data is not “long” enough to uniquely determine S(0); that is, almost the same quality of data fit can be achieved for different numbers S(0) (we have tried 50,000, 100,000 and 200,000) leading to different numbers of “stabilized” cumulative infected cases and infection periods. Taking this factor into account, we do not consider this approach, however we note that it might be quite possible soon with the availability of more data points.

3. In this paper we adopt another approach by neglecting the compartment S completely and leaving just the compartment I. The force of infection F(S, I) in this case is the main factor to be determined. We take this function in the form

F(S, I) = β𝒟I           (5)

where 𝒟 is the population density of a particular country. In a more general setting, one would involve functions nonlinear in I (like F(S, I) = β𝒟I^ξ with ξ ≤ 1). However, since the infectious population I is a very small portion of the total population, function F can be assumed linear at least in early stages of epidemics. In this case equation (4) can be represented in the form

I(t + 1) = (1 – μ)^τ₁ β𝒟I(t – τ₁) + (1– μ – α – γ)I(t).           (6)

The major drawback of this model is that I may growth infinitely if the reproduction number is greater than 1; in this model there is no variable/parameter (like S(t) in SIR) that could force I to decrease. On the other hand we believe that it can better describe the behavior of an infected population in “small” time intervals and provide more accurate reproduction numbers.

Active infectious population. Now we discuss the infectious population and equation (6) in more detail. We call “active infectious populations” at time t the infected population that are infectious at that time but are not hospitalized yet. Denote by I_a(t) the number of active infectious populations at time t. We will rewrite equation (6) in terms of I_a.

Denote by τ₂ the average infectious period; that is, time from onset (τ₁) to hospitalization. Then, an infected person is assumed to be active infectious during the period [τ₁, τ₁ + τ₂]. Since τ₂ is relatively small, we can assume that none is recovering during that period. This means that the rate of recovery γ in (6) is no longer needed for I_a(t).

Thus, we transform equation (6) by taking into account the time delay τ₂. Accordingly, the equation for I_a(t) can be represented in the form

Here ω(0) = 0 and ω(i), i ≥ 1, is a gamma cumulative distribution function for onset-to-death that well describes the current Ebola virus in West Africa⁷. We note that in this equation, for each i ≥ 1, the fraction (1–αω(i)) is applied to the remaining infectious (1–μ)ⁱ β𝒟I_a(t–τ₁–i); that is, the death rate in (7) is slightly different from (6) (indeed, both μ and ω(i) are quite small and this leads to 1–μ–αω(i) ≈ (1–μ)(1–αω(i))).

Cumulative number of infected cases. The first term (1 – μ)^τ₁ β𝒟 I_a(t – τ₁) in (7) describes the number of new cases at time t. The cumulative number of infectious cases at (t+1) will be denoted by C(t +1). It can be calculated as

Cumulative number of deaths. To calculate the cumulative number of deaths at time t, we consider all infectious cases (hospitalized or not) in the interval [t–τ₁, t–n] where n is a sufficiently large number. In particular we assume that death may occur after the onset. As mentioned above, the distribution of death is described by a gamma distribution function ω with its p.d.f - ω_p. Then, the cumulative number of deaths due to disease can be calculated as