Be-CoDiS: A mathematical model to predict the risk of human diseases spread between countries. Validation and application to the 2014-15 Ebola Virus Disease epidemic.

Be-CoDiS: A mathematical model to predict the risk of human diseases spread between countries. Validation and application to the 2014-15 Ebola Virus Disease epidemic.

Benjamin Ivorra, Diène Ngom and Ángel M. Ramos
 
MOMAT research group, IMI-Institute and
Applied Mathematics Department.
Complutense University of Madrid.
Plaza de Ciencias, 3, 28040, Madrid, Spain.
 
Département de Mathématiques
UFR des Sciences et Technologies
Université Assane Seck de Ziguinchor &
Laboratoire d’Analyse Numérique et d’ Informatique,
Université Gaston Berger de Saint Louis, Sénégal.
 
Corresponding author: Tel: +34 91 394 44 15;
E-mail: ivorra@mat.ucm.es
Version 5.0 - Date of publication: July 13, 2019
Abstract

Ebola virus disease is a lethal human and primate disease that currently requires a particular attention from the international health authorities due to important outbreaks in some Western African countries and isolated cases in the United Kingdom, the USA and Spain. Regarding the emergency of this situation, there is a need of development of decision tools, such as mathematical models, to assist the authorities to focus their efforts in important factors to eradicate Ebola. In this work, we propose a novel deterministic spatial-temporal model, called Be-CoDiS (Between-Countries Disease Spread), to study the evolution of human diseases within and between countries. The main interesting characteristics of Be-CoDiS are the consideration of the movement of people between countries, the control measure effects and the use of time dependent coefficients adapted to each country. First, we focus on the mathematical formulation of each component of the model and explain how its parameters and inputs are obtained. Then, in order to validate our approach, we consider two numerical experiments regarding the 2014-15 Ebola epidemic. The first one studies the ability of the model in predicting the EVD evolution between countries starting from the index cases in Guinea in December 2013. The second one consists of forecasting the evolution of the epidemic by using some recent data. The results obtained with Be-CoDiS are compared to real data and other models outputs found in the literature. Finally, a brief parameter sensitivity analysis is done. A free Matlab version of Be-CoDiS is available at: http://www.mat.ucm.es/momat/software.htm

Keywords: Epidemiological modelling; Ebola Virus Disease; Deterministic models; Compartmental models

1 Introduction

Modeling and simulation are important decision tools that can be used to control or eradicate human and animal diseases (Anderson, 1979; Thieme, 2003). Each disease presents its own characteristics and, thus, most of them require a well-adapted simulation model in order to tackle real situations (Bowman et al, 2005; Brauer and Castillo-Chávez, 2001).

In this work, we present a first version of a new deterministic spatial-temporal epidemiological model, called Be-CoDiS (Between-Countries Disease Spread), to simulate the spread of human diseases in a considered area. This model is inspired from a previous one, called Be-FAST (Between Farm Animal Spatial Transmission), which focuses on the spread of animal diseases between and within farms. The major original ideas introduced by Be-FAST were the following: (i) the study of both within and between farms spread, (ii) the use of real database and (iii) dynamic coefficients calibrated in time according to farms characteristics (e.g., size, type of production, etc.). This model was deeply detailed in Ivorra et al (2014); Martínez-López et al (2011) and validated on various cases as, for example: Classical Swine Fever in Spain and Bulgaria (Martínez-López et al, 2012, 2013), and Foot-and-Mouth disease in Peru (Martínez-López et al, 2014). Be-CoDiS is based on the combination of a deterministic Individual-Based model (with the countries playing the role of individuals) (DeAngelis and Gross, 1992), simulating the between-country interactions (here, movement of people between countries) and disease spread, with a deterministic compartmental model (Brauer and Castillo-Chávez, 2001) (a system of ordinary differential equations), simulating the within-country disease spread. We observe that the coefficients of the model are calibrated dynamically according to the country indicators (e.g., economic situation, climatic conditions, etc.). At the end of a simulation, Be-CoDiS returns outputs referring to outbreaks characteristics (for instance, the basic reproductive ratio, the epidemic magnitude, the risk of disease introduction or diffusion per country, the probability of having at least one infected per time unit, etc.). The main characteristics of our approach are the consideration of movement of people between countries and control measure effects and the use of time dynamic coefficients fitted to each country.

This work has two main goals. The first one, is to give a full and detailed mathematical formulation of Be-CoDiS and to explain how the model parameters are obtained in order to provide a transparent and understandable model for users. The second one is to validate our model by applying it to the current case of the Ebola Virus Disease (EVD) (Emond et al, 1977; Peters and Peters, 1999; W.H.O., 2014a).

EVD is a lethal human and primates disease caused by the Ebola-virus (family of the Filoviridae) that causes important clinical signs, such as haemorrhages, fever or muscle pain. The fatality percentage (i.e., the percentage of infected persons who do not survive the disease) is evaluated to be within 25% and 90%, due to hypovolemic shock and multisystem organ failure, depending on the sanitary condition of the patient and the medical treatment. This virus was first identified in Sudan and Zaire in 1976 (see Emond et al (1977)). Various important outbreaks occurred in 1995, 2003 and 2007 in the Democratic Republic of Congo (315, 143 and 103 persons infected by EVD, respectively), in 2000 and 2007 in Uganda (425 and 149 persons infected by EVD, respectively), see Peters and Peters (1999); Chowell et al (2004). The 2014-15 outbreak started in December 2013 in Guinea and spread to Liberia and Sierra Leone. In March 2014, the international community was aware of the gravity of the situation in those three countries. The situation on April 24, 2015 (the date used to run our numerical experiments) was a total of 26307 persons infected by EVD in Guinea, Liberia, Sierra Leone and Nigeria (see W.H.O. (2014a, b)). Moreover, 15 isolated cases were detected in Mali, Senegal, Spain, the United Kingdom and the USA. Furthermore, in Spain, the United Kingdom and the USA, the first contagions between people outside Africa were observed. The observed mean fatality percentage for this particular hazard decreased from 72.8% (on March, 2014) to 47.5% (on April, 2015), see Leroy et al (2004). From an epidemiological point of view, the EVD can be transmitted between natural reservoirs (for instance, bats) and humans due to the contact with animal carcass (Leroy et al, 2004). The most common way of EVD human transmission is due to contacts with blood or bodily fluids from an infected person (including dead persons).

Starting from this particular context, we study the behaviour of our model in predicting the possible spread of EVD worldwide. In order to validate our approach, we first consider a numerical experiment starting from the index cases in Guinea in December 2013 and check the ability of the model in spreading the EVD to other countries. Then, we perform a second experiment which aims to forecast the possible evolution of EVD by considering recent data. Moreover, we also simulate the possible evolution of EVD in the countries with the highest risks of EVD introduction due to the movement of people. The results obtained by Be-CoDiS at the end of those experiments are compared to historical data and with other studies found in the literature C.D.C. (2014); Fisman et al (2014); Gomes et al (2014); Meltzer et al (2014); W.H.O. (2014a, b); W.H.O. Response Team (2014) regarding the same 2014-15 EVD outbreaks. The work in Gomes et al (2014) is based on a time spatial model with stochastic flux between areas (considering a database based on airport traffic instead of, as done here, calibrated values of migratory flux) but with a different point of view regarding control measures (at the beginning of the simulation, if control measures are applied the authors set the disease transmission rate to a value lower than in the case without control measures) and model coefficients (considered as constant in time). We point out that most of the parameters used by our model are calibrated for African countries and few data are available about the behaviour of EVD in other countries. Thus, some empirical hypothesis are needed and a sensitivity analysis of Be-CoDiS regarding those parameter is done. Finally, we also highlight the current limitations of the model and a way to improve it in future works.

This work is organized as follows. In Section 2, we describe the epidemiological behavior of EVD and give a detailed presentation of our model. In particular, we focus on its mathematical formulation and explain how its parameters are obtained. In Section 3, we focus on the validation of Be-CoDiS by simulating possible evolutions of the EVD spread worldwide and by comparing the obtained results with observed data and other studies. Finally, we study the behavior of our model regarding changes in the whole set of values of its parameters.

2 Mathematical formulation of the model

In Section 2.1 we detail the epidemiological characteristics of EVD that are taken into account in our model. Then, in Sections 2.2, 2.3 and 2.4, we describe in detail the Be-CoDiS model by presenting its general structure, i.e., the considered within and between countries disease spread sub-models. Finally, in Section 2.5, we present the outputs used to analyse the results of the numerical simulations presented in Section 3. The main notations used in this work are summarized in Table LABEL:tnota.

Remark 1

Be-CoDiS is designed to be able to study the spread of any human disease worldwide. Here, some particular details of the model are related to the specific EVD case but it can be adapted to other disease. For instance, the classification of compartments in the SEIHRDB model can be changed to study other cases.

Notation Value Description
[0,+) Maximum number of simulation days (day)
1 Time discretisation step size (day)
[0.1657,0.2671] Disease contact rate of a person
in state in country (day)
Disease contact rate of a person
in state in country (day)
Disease contact rate of a person
in state in country (day)
0.0877 Transition rate of a person in state (day)
[0.2,0.5] Transition rate of a person in state (day)
in country at time ,
[0.14,0.2] Transition rate of a person in state
to state (day) in country at time ,
[0.13,0.24] Transition rate of a person in state
to state (day) in country at time ,
12.9 the period of convalescence (day)
[0.012,0.023] Natural mortality rate in country (day)
[0.22,1.37] Natural natality rate in country (day)
[0.25,0.728] Disease fatality percentage in country at time
[0.001,0.281] Efficiency of the control measures
in country (day)
[0,+) First day of application of
control measures in country (day)
[0,1] Control measure efficiency (%) in country at
time applied to persons in state or
[0,1] Control measure efficiency (%) applied to persons
moving from country to country at time
Daily rate of the movement of people
from country to country (day)
0.53 Proportion of the that can be reduced
due to the application of control measures
176 Number of countries
[2.5,1.4] Number of persons in country at time
[2.5,1.4] Number of persons in state , , , , , ,
in country at time
Table 1: Summary of the main notations used in this work to describe Be-CoDiS. A brief description (Description) and the range of the considered values (Value) are also reported. The reference are given in the corresponding section of the text.

2.1 Epidemiological characteristics of Ebola Virus Disease

When a person is not infected by EVD, it is categorized in the Susceptible state (denoted by ). If a person is infected, then they passes successively through the following states (see Legrand et al (2007); Meltzer et al (2014); Peters and Peters (1999); W.H.O. Response Team (2014)):

  • Infected (denoted by ): The person is infected by EVD but they cannot infect other people and has no visible clinical signs (i.e., fever, hemorrages, etc.). The mean duration of a person in this state is 11.4 days (range of [2,21] days) and is called incubation period. Then, the person passes to the infectious state.

  • Infectious (denoted by ): The person can infect other people and start developing clinical signs. The mean duration of a person in this state is called infectious period. In September, 2014, the mean infectious period was 5 days (range of [0,10] days, according to W.H.O. Response Team (2014)). However, due to the control measures applied to fight the EVD epidemic presented below, in March, 2015, the mean infectious period is estimated to be between [2.6,3] days for affected African countries (see, W.H.O. (2014b)). Furthermore it took only 2 days to hospitalize an infected patient in the United Kingdom. After this period, infectious persons are taken in charge by sanitary authorities and we classify them as Hospitalized.

  • Hospitalized (denoted by ): The person is hospitalized and can still infect other people, but with a lower probability. The mean duration in this sate is 4.5 days (see Legrand et al (2007)). Actually, it has been observed in practice that those patients can still infect other people with a probability 25 times lower than infectious people (see W.H.O. (2014b)). On the one hand, after a mean duration of 4.2 days (range [1,11] days), a percentage of the hospitalized persons, between 25% and 72.8% (depending on the sanitary services of the country), die due to the EVD clinical signs and pass to the Dead state. On the other hand, after a mean duration of 5 days, the persons who have survived to EVD pass to the Recovered state. We point out that state does not include hospitalized persons which cannot infect other people any more. This last category of persons is included in the recovered state explained below. The mean of the total number of days that a EVD patient stays in hospital (from hospitalization to hospital discharge) is estimated to be 11.8 days (range [7.7,17.9] days).

  • Dead (denoted by ): The person has not survived to the EVD. It has been observed in previous Ebola epidemics, that cadavers of infected persons can infect other people until they are buried. The probability to be infected by this kind of contact is the same as that of being infected by contact with infectious persons. Indeed, the daily number of contacts of a cadaver with persons is lower than that of an infectious person but, the risk of EVD transmission due to contacts with cadavers is larger than the risk due to contacts with infectious persons (see Legrand et al (2007)). After a mean period of 2 days (this could be different for some countries) the body is buried.

  • Buried (denoted by ): The person is dead because of EVD. Its cadaver is buried and they can no longer infect other people.

  • Recovered (denoted by ): The person has survived the EVD and is no longer infectious. They develop a natural immunity to EVD. Since it has never been observed a person who has recovered from Ebola and contracted the disease again during the period of time of the same epidemic, it is assumed that they cannot be infected again by Ebola.

Once an EVD infected person is hospitalized, the authorities apply various control measures in order to control the EVD spread (see Fisman et al (2014); Hernandez-Ceron et al (2013); W.H.O. (2014a, b)):

  • Isolation: Infected people are isolated from contact with other people. Only sanitary professionals are in contact with them. However, contamination of those professionals also occur (see Fisman et al (2014)). Isolated persons receive an adequate medical treatment that reduces the EVD fatality rate.

  • Quarantine: Movement of people in the area of origin of an infected person is restricted an controlled (e.g., quick sanitary check-points at the airports) to avoid that possible infected persons spread the disease.

  • Tracing: The objective of tracing is to identify potential infectious contacts which may have infected a person or spread EVD to other people.

  • Increase of sanitary resources: The number of operational beds and sanitary personal available to detect and treat affected persons is increased, producing a decrease in the infectious period. This greatly the time needed to The funerals of infected cadavers are controlled by sanitary personal in order to reduce the contacts between the dead bodies and susceptible persons.

Remark 2

Data given above were calibrated for the cases of African countries, the natural reservoir of EVD. However, due to the spread of this disease out of Africa, new studies should be performed to analyze the behavior of EVD in other sanitary, population and climatic conditions. Currently, very few studies are available. One of them is about the survival of the Ebola virus (EV) according to changes in temperature (see Piercy et al (2010)). It has been found that the lower is the temperature the greater is the survival period of the EV) outside the host. Thus, in this work some empirical hypothesis, which seems to be reasonable, have been done. We have assumed that the transmission parameter of EV decreases when the temperature or the sanitary expenses of a country increase and increases when the people density of a country increases.

2.2 General description

The Be-CoDiS model is used to evaluate the spread of a human disease within and between countries during a fixed time interval.

At the beginning of the simulation, the model parameters are set by the user (for instance, the values considered for EVD are described in Section 3.1). At the initial time (), only susceptible people live in the countries that are free of the disease, whereas the number of persons in states , , , , , and of the infected countries are set to their corresponding values. Then, during the time interval , with being the maximum number of simulation days, the within-country and between-country daily spread procedures (described in Section 2.3) are applied. If at the end of a simulation day all the people in all the considered countries are in the susceptible state, the simulation is stopped. Else, the simulation is stopped when . Furthermore, the control measures are also implemented and they can be activated or deactivated, when starting the model, in order to quantify their effectiveness to reduce the magnitude and duration of a EVD epidemic.

A diagram summarizing the main structure of our model is presented in Figure 1.

Remark 3

The choice of using a deterministic model instead of a stochastic one is done as a first approach, since such kind of models presents some advantages, such as: a low computational complexity allowing a better calibration of the model parameters or the possibility of using the theory of ordinary differential equations for well analyzing and interpreting the model. Furthermore, according to Diekmann et al (2012), deterministic models should be the first tool to be used when modelling a new problem with few data (here, few data of the spread of Ebola are available for non-African countries). The authors of that work also note that the stochastic models are not suitable when it is difficult or impossible to determine the distribution probability, are difficult to analyze and require more data for the calibration of the model.

Figure 1: Diagram summarizing the Be-CoDiS model.

2.3 Within-Country disease spread

The dynamic disease spread within a particular contaminated country is modeled by using a deterministic compartmental model (see, for instance, Brauer and Castillo-Chávez (2001)).

We assume that the people in a country are characterized to be in one of those states, described in Section 2.1: Susceptible (), Infected (), Infectious (), Hospitalized (), Recovered (), Dead () or Buried (). For the sake of simplicity we assume that, at each time, the population inside a country is homogeneously distributed (this can be improved by dividing some countries into a set of smaller regions with similar characteristics). Thus, the spatial distribution of the epidemic inside a country can be omitted. We also assume that new births are susceptible persons. In Section 2.3, we do not consider interaction between countries.

Under those assumptions, the evolution of , , , , , and , denoting the number of susceptible, infected, infectious, hospitalized, recovered, dead and buried persons in country at time , respectively, is modeled by the following system of ordinary differential equations

(1)

where:

  • ,

  • is the number of countries,

  • is the number of persons (alive and also died or buried because of EVD) in country at time ,

  • is the natality rate (day) in country : the number of births per day and per capita,

  • is the mortality rate (day) in country : the number of deaths per day and per capita (or, equivalently, the inverse of the mean life expectancy (day) of a person),

  • if , if , and 0 elsewhere, with (small tolerance parameter). This function is a filter used to avoid artificial spread of the epidemic due to negligible values of ,

  • is the disease fatality percentage in country at time : the percentage of persons who do not survive the disease,

  • is the disease effective contact rate (day) of a person in state in country : the mean number of effective contacts (i.e., contacts sufficient to transmit the disease) of a person in state per day before applying control measures,

  • is the disease effective contact rate (day) of a person in state in country ,

  • is the disease effective contact rate (day) of a person in state in country ,

  • , , , , denote the transition rate (day) from a person in state , , , or to state , , , or , respectively: the number of persons per day and per capita passing from one state to the other (or, equivalently, the inverse of the mean duration of one of those persons in state , , , , or , respectively). We note that , and are time and country dependent, since, due to the applied control measures in country , their value could vary in time,

  • , , (%) are functions representing the efficiency of the control measures applied to non-hospitalized persons, hospitalized persons and infected cadavers respectively, in country at time to eradicate the outbreaks. Focusing on the application of the control measures, which in the EVD consists in isolating persons/areas of risks and in improving sanitary conditions of funerals, we multiply the disease contact rates (i.e., , and ) by decreasing functions simulating the reduction of the number of effective contacts as the control measures efficiency is improved. Here, we have considered the functions (see Lekone and Finkenstädt (2006)):

    (2)

    where (day) simulates the efficiency of the control measures (greater value implies lower value of disease contact rates) and (day) denotes the first day of application of those control measures.

System (1) is completed with initial data , , , , , and given in , for =1,.., .

We observe that all parameters of system (1) should be adapted to the considered disease and countries. Generally, they are calibrated considering real data as explained in Section 3 for the EVD case.

Remark 4

Numerical experiments presented in Section 3.2 seem to show that, despite its apparent simplicity, using system (1) to simulate the within-country disease spread gives reasonable results. However, a next step could be to model the wihtin-country spread by considering, for instance, an individual based model that simulates the interactions between communities (e.g., cities, villages, etc.), as done in the case of herds and livestock diseases in Martínez-López et al (2014). However, to do so, we first need to collect data about those communities (i.e., size, location, type, etc.) and their interacting network (i.e., movement of people, etc.). Obtaining those data at the worldwide level seems to be quite difficult and may require a large effort of data-mining. Furthermore, the computational time needed to solve such a model would be much higher than that of the current version of Be-CoDiS and, thus, the estimation of the model parameters would be complicated. However, some of the interests in developing this kind of complex model are, for instance, their ability to estimate the local efficiency of the control measures and to assess the geographical allocation of the control measures resources in order to contain efficiently the epidemic spread.

2.4 Between-Countries disease spread

The disease spread between countries is modeled by using a spatial deterministic Individual-Based model (see DeAngelis and Gross (1992)). Countries are classified in one of the following states: free of disease () or with outbreaks (). We assume that at time country is in state if (i.e., there exist at least one infected person in this country), else it is in state .

In this work we consider that the flow of people between countries and at time (i.e., persons traveling per day from to at time ), is the only way to introduce the disease from country , in state , to country . To do so, we consider the matrix , where is the rate of transfer (day) of persons from country to country , which is expressed in % of population in per unit of time. Furthermore, we assume that only persons in the and sates can travel, as other categories are not in condition to perform trips due to the importance of clinical signs or to quarantine. Moreover, due to control measures in both and countries, we assume that those rates can vary in time and are multiplied by a function denoted by .

Thus, we consider the following modified version of system (1):

(3)

System (3) is completed with initial data , , , , , and given in ; for =1,.., .

This full model (3), which is summarized in Figure 1, is called Be-CoDiS.

Remark 5

The explicit solution of (3) and the corresponding initial values are not (in general) available. In order to get an approximation of the solution one can use a suitable numerical solver. Here, we use the explicit Euler scheme with a time step of 1 day.

2.5 Outputs of the model

Here, we present the outputs used to analyse the results of the simulations performed in Section 3. In particular, considering a time interval , for each country we compute the following values:

  • : the cumulative number of EVD cases in country at day , which can be computed as:

  • : cumulative number of deaths (due to EVD) in country at day , which can be computed as:

  • R and CR: the basic reproductive ratio of the model and country , respectively. It is defined as the number of cases one infected person generates on average over the course of its infectious period, in an otherwise uninfected population. In our model most of the parameters are time and country dependent, thus in order to approximate a value of R we propose the following methodology. First, considering the approach proposed in Gomes et al (2014); Legrand et al (2007), we compute , an estimation of the the R value obtained when considering system (1) with the parameters value of country at time , given by

    Then, according to the number of infected persons in each country at time , we compute the mean basic reproductive ratio of the model at time , denoted by , as

    Finally, we compute

    (4)
  • RS: the initial risk of country to spread EVD to other countries, given by:

    RS (personsday) is the daily amount of infected persons who leave country at time .

  • TRS: the total risk of country to spread EVD to other countries, considering the time interval , computed as:

    TRS (persons) is the number of infected persons send to other countries during the considered time interval.

  • RI: the initial risk of EVD introduction into country from other countries, given by:

    RI (personsday) is the daily amount of infected persons entering country at time .

  • TRI(: the total risk of EVD introduction into country from other countries considering the time interval , computed as:

    TRI (persons) is the number of EVD infected persons received from other countries during the time interval .

  • MNH(): the maximum number of hospitalized persons at the same time at country during the time interval . It is computed as:

    We remind (see Section 2.1) that is the period of convalescence (i.e., the time a person is still hospitalized after surviving EVD). This number can help to estimate and plan the number of clinical beds needed to treat all the EVD cases.

  • Emf: the percentage of the population leaving the country each day (day), which can be computed as:

Remark 6

We note that and are some of the main indicators reported in W.H.O. (2014b), since March 2014, in order to give an estimation of the magnitude of the epidemic.

3 Application to the 2014-15 Ebola case

We are now interested in validating our approach by considering the Ebola epidemic currently occurring worldwide. The advantage of this case is that some real and simulated data are now available and, thus, we are able to compare our model outputs with the information available in the following literature Fisman et al (2014); Astacio et al (2014); W.H.O. (2014a).

To this aim, we first explain in Section 3.1 how to estimate the model parameters for the EVD. Then, we present the results and discuss them in Section 3.2. Finally, in Section 3.3, we carry out a brief sensitivity analysis regarding the parameters values.

3.1 Be-CoDiS parameters estimation for EVD

Some of the parameters used in the simulations presented in Section 3.2, have been found in the literature Fisman et al (2014); Astacio et al (2014); Chowell et al (2004); Gomes et al (2014) and in the daily reports on the Ebola evolution available online (see C.D.C. (2014); W.H.O. (2014b)). Despite the effort to use the maximum amount of robust parameters as possible, due to lack of information of the behavior of Ebola out of Africa, some of them have been estimated using empirical assumptions. This part should be clearly improved as soon as missing information is available.

We now detail each kind of parameter by its category.

3.1.1 Country indicators

We use the following data regarding country :

  • : Mean temperature (ºC) at day .

  • : Population density (persons/km).

  • : Total number of persons alive and also died or buried because of EVD at time (see W.H.O. (2014b)).

  • : Gross National Income per year per capita (US.person.ye-ar). We remind that the Gross National Income is an indicator of the economy of the country: the total domestic and foreign output claimed by residents of a country.

  • : The mean Health expenditure per year per capita (US.per-son.year). This is an economical indicator of the sanitary system of a country given by the amount of money inverted by public and private institutions (national or international) in the sanitary system of the country.

  • : The mean life expectancy (days).

  • : See Section 2.3.

All those data have been freely obtained for year 2013 from the following World Data Bank website: http://data.worldbank.org.

3.1.2 Initial conditions

We have considered the initial conditions for our system (3) corresponding to the state of the EVD epidemic at several dates reported in W.H.O. (2014b).

In W.H.O. (2014b), the cumulative numbers of reported cases (i.e., persons who have ever been in state , as stated in W.H.O. Response Team (2014)) and deaths (i.e., a person in state or ) in each country at date , denoted by NRC and NRD, are given for dates starting on March, 23 2014. However, those report are not published daily. Thus, in order to complete the missing information we use a cubic hermite interpolation method assuming that on November 1, 2013 all countries were free of the disease. Thus, taking into account the characteristics of the EVD presented in Section 2.1, we estimate the amount of persons in states , , , , and at as follows:

  • Since the main duration in state is 4.5 days, we compute by considering the number of reported cases that are alive at time minus the number of reported cases that are alive 4.5 days before

  • Since the main duration in state is 2 days, we compute as

  • Since the main duration in states , and is 11.4, 5 and 4.5 days, respectively, we consider

    where is the total number of persons who have ever been in state any of the last 4.5 days.

  • The number of recovered persons is given by

  • The number of buried cadaver is

  • The number of susceptible persons is

All these numbers are rounded to the nearest integer.

3.1.3 Rate of movement of people between countries

Regarding the dynamic of the current 2014-15 Ebola epidemic and the way the disease has spread from Guinea to other countries, it seems that the short-term visits (such as the visit to relatives) are the main causes of EVD diffusion between countries (Chowell and Nishiura, 2014). However, such information is quite difficult to obtain at a worldwide level, especially when considering African countries.

Here, to approximate the pattern of those short term visits, we have considered the data regarding the 2005-2010 human migratory fluxes between countries obtained from Guy and Nikola (2014) and freely available at the URL: http://www.global-migration.info Following the analysis provided in O.N.S. (2013), that compare short term immigration (a stay round 3 months) and the census in England, we assume that the pattern of the movement of people between countries of this migration database and the short term visits exhibit a similar general behavior (i.e., similar visited countries). However, it is expected that short term visits occur much more often than long-term immigration. Thus, considering both assumptions, we compute the percentage of persons in country moving to country per day as

where represents the number of persons moving from country to country from 2005 to 2010 and is a parameter used to increase the value of .

In addition to , the parameter , presented in Section 2.4, also play a role in the control of the movement of people as it does not allow a country to spread EVD to other countries until . In this work, we have estimated, by considering a trial and error process performed in the experiment presented in Section 3.2.2, that and produce a pattern of EVD spread between countries close to the one observed during the 2014-15 epidemic. This numerical result seems to indicate that the approach considered here regarding the movement of people is reasonable.

Of course, if data about short term visits are available they should replace the values of matrix used in this work.

3.1.4 EVD characteristics

The following parameters are assumed to be well studied due to several data sets available about the current 2014-15 Ebola outbreak. Using data from Sections 2.1 and 3.1.1, we estimate the following parameters for our model:

  • (day)

  • , where we denote by W the disease fatality percentage of country when no control measure is applied; is the minimum disease fatality percentage; is the maximum disease fatality percentage; and denotes the proportion of the fatality percentage that can be reduced due to the application of control measures. For the EVD, we consider , and . This formula and the value of were determined empirically during the experiments presented in Section 3.2 so that the cumulative number of deaths returned by our model and the real observations reported in W.H.O. (2014b) were similar. To determine the expression of , we have taken into account that, for this particular epidemic, the EVD fatality percentage oscillate between [25,72.8]%, depending on the quality of the sanitary service (see Meltzer et al (2014)). Moreover, according to W.H.O. (2014b), the maximum fatality rate has decreased with the application of the control measures. Thus, we have modelled this effect by multiplying the maximum fatality rate by .

  • (day) and (day)

  • 13 (days) is an upper bound estimation of the convalescence period computed as the maximum number of days that a persons stays hospitalized minus the duration of a person in the state ,

  • (day), (day), (day), where , and denote the mean duration in days of a person from state I to H, from state H to R and from state H to D, respectively, without the application of control measures; represents the decrease of the duration due to the application of control measures in country at time ; and is the maximum number of days that can be decreased due to the control measures. Here, according to Section 2.1, we consider , , and (assuming that in cases of strong control measures, as in the United Kingdom, can be reduced by 3 days, see W.H.O. (2014b)).

  • : There exists several works on the computation of the EVD effective contact rate considering various SIR model (see Astacio et al (2014); Chowell et al (2004)). However, the value of this rate depend on the epidemic characteristics (country, year, etc.). Furthermore, our model includes novel characteristics regarding those articles, as it includes movement between countries, hospitalized people and control measures. Thus, we have computed our own rates by using a regression method considering three particular sets of data associated with the evolution of the EVD epidemic in Guinea, Liberia and Sierra Leone (see W.H.O. (2014b); C.D.C. (2014)).

    • In Guinea, the country of origin of the EVD epidemic, the index case was identified as a young boy who died on December 6, 2013 and infected 3 persons of its family. On July 14, 2014, the estimated date for which the effects of the application of the control measures by national and international authorities start to slow down the EVD spread dynamic, a total cumulative number of 425 cases and 319 dead persons were reported (see W.H.O. (2014b)). After this date, the international help started to affect the initial EVD effective contact rate in Guinea, denoted by . Thus, we fit those data with the solution given by system (1). To this end, system (1) was started at (corresponding to December 6, 2013) with 3 persons in state in Guinea, 1 person in state and all other persons being free of disease. The model was run with days (corresponding to July 14, 2014 as the final date). In this particular simulation, we did not consider control measures (i.e., for any , ). All other parameters were set to the values introduced previously. Considering a particular value , at the end of this simulation we computed the model error Err. We minimized Err by considering a dichotomy algorithm starting from =0.117 (day) (see Astacio et al (2014)) and found an optimal value of =0.1944 (day).

    • In Sierra Leone, 8 cases and 2 deaths were reported on March 31, 2014. The effects of control measures are estimated to start on September 9, 2014,when the cumulative number of reported cases was 2792 with 1501 deaths. We used the same fitting technique as in the case of Guinea, with system (1), without control measures, starting at (corresponding to March 31, 2014) with 3, 1, 1, 1, 5 and 1 persons in states , , , , and , respectively (considering the estimation method presented in Section 3.1.2) and all other persons being free of disease. The system was run with days (i.e., final date on September 9, 2014). We found =0.2605 (day).

    • In Liberia, 16 cases and 5 deaths were reported on May 27, 2014. On June 30, 2014, before the start of effects of the control measures, 302 cases and 139 deaths were observed. System (1), without control measures, was started at (i.e., corresponding to May 27, 2014) with 5, 2, 2, 1, 9 and 4 persons in the states , , , , and , respectively (see Section 3.1.2). This system was run during days. Applying the same technique as for Guinea and Sierra Leone, we found =0.2649 (day).

    Taking into account those three rates, since the rate of other countries (especially the non African ones) remains unknown (due to the lack of data), we have performed an empirical non linear regression to estimate . To this aim, is assumed to be a non-decreasing function , where (kmUS$personsyear) is a balance parameter which determines the importance of on the value of in comparison to . Indeed, the variable is chosen because of the following reasons: 1) we assume that the higher the population of a country is, the higher the probability of contagion is and the higher is; 2) the higher the economy level of a country is, the higher its education level is, the lower the EVD risk habits of persons are (for instance, touching cadavers during funerals, see Fonkwo (2008); Hewlett and Hewlett (2007)) and the lower is. In addition, we propose to use a function of the form

    where (day ), (non-dimensional) and (day) . We found, by considering the nonlinear regression method ’nlinfit’ implemented in Matlab using the points

    • ,

    • Sierra Leone,

    • Liberia,

    that , , and .

    Remark 7

    If needed, can be also considered as time dependent (see Forgoston and Schwartz (2013)). For instance, as said in Section 2.1, it has been observed that Ebola Virus survives better outside the host for lower temperatures (see Piercy et al (2010)). Thus, it could be interesting to introduce a slight dependence of on the temperature of the country . For instance, we could consider:

    (5)

    where TMP (ºC) is a reference temperature; and (%) represent the maximum percent variation of the value . However, as no data are available in literature to estimate a suitable value of , the effect of the temperature is neglected in our model.

  • (day), since the probability of being infected by contact with persons in state is 25 times lower than the probability of being infected by contact with persons in state , as explained in Section 2.1,

  • (day), since the probability of being infected by contact with persons in state is the same as the probability of being infected by contact with persons in state (see Section 2.1).

3.1.5 Control measures

Here, we estimate the parameters used in Equation 2:

  • : It is the first day such that for all countries except for Guinea, Liberia, Mali, Nigeria and Sierra Leone. For these countries, intensive control measures were not applied right after the apparition of the first person in state , but some time later (as reported in W.H.O. (2014b)). Thus, in the simulations presented in Section 3.2, we considered for these countries a reported delay between the detection of the first EVD case and the application of the control measures.

  • : In order to fit , we used data from Guinea, Sierra Leone and Liberia. Moreover, for all the numerical experiments considered in this work, we considered .

    • In Guinea, the effect of the control measures started on July 14, 2014. On April 15, 2015, the number of reported cases in Guinea was 3568. Again, we fit those data with system (1) starting at (corresponding to December 6, 2013) with 3 persons in state and 1 person in state in Guinea and all other persons being free of the disease. The model was run with =496 days (corresponding to April 15, 2015). In this simulation, the control measures were applied after days (i.e., (Guinea)=220 days). Considering a particular value , we computed the model error Err, as defined previously. We minimized Err by considering a dichotomy algorithm starting from =0.001 and found an optimal value of =0.00125 (day).

    • In Sierra Leone, the effects of control measures started on September 9, 2014. On April 15, 2015, the number of reported cases in Sierra Leone was 12294. We used the same fitting method as in Guinea and started system (1) with the same conditions as those used for computing . The system was started with =389 days and control measures were applied at day (Sierra Leone)=228 (i.e., corresponding to September 9, 2014). We found =0.00227.

    • In Liberia, the effects of control measures started on June 30, 2014. On April 15, 2015, the number of reported cases was 10241. We used the same fitting method as the one used for Guinea and started system (1) with the same conditions as those used for computing . This system was run with =333 days and control measures were applied at day =34. We found =0.00270.

    Taking into account those three values, we perform a regression method, similar to the one introduced in Section 3.1.4, estimating . To this end, is assumed to be a non-decreasing function , where (personsyear kmUS$) is a balance parameter which determines the importance of