A Mathematical Model with Quarantine States for the Dynamics of Ebola Virus Disease in Human Populations.

A deterministic ordinary differential equation model for the dynamics and spread of Ebola Virus Disease is derived and studied. The model contains quarantine and nonquarantine states and can be used to evaluate transmission both in treatment centres and in the community. Possible sources of exposure to infection, including cadavers of Ebola Virus victims, are included in the model derivation and analysis. Our model's results show that there exists a threshold parameter, R 0, with the property that when its value is above unity, an endemic equilibrium exists whose value and size are determined by the size of this threshold parameter, and when its value is less than unity, the infection does not spread into the community. The equilibrium state, when it exists, is locally and asymptotically stable with oscillatory returns to the equilibrium point. The basic reproduction number, R 0, is shown to be strongly dependent on the initial response of the emergency services to suspected cases of Ebola infection. When intervention measures such as quarantining are instituted fully at the beginning, the value of the reproduction number reduces and any further infections can only occur at the treatment centres. Effective control measures, to reduce R 0 to values below unity, are discussed.

1. Introduction and Background

The world has been riveted by the 2014 outbreak of the Ebola Virus Disease (EVD) that affected parts of West Africa with Guinea, Liberia, and Sierra Leone being the most hard hit areas. Isolated cases of the disease did spread by land to Senegal and Mali (localized transmission) and by air to Nigeria. Some Ebola infected humans were transported to the US (except the one case that traveled to Texas and later on died) and other European countries for treatment. An isolated case occurred in Spain, another in Italy (a returning volunteer health care worker), and a few cases in the US and the UK [1-3]. Though dubbed the West African Ebola outbreak, the movement of patients and humans between countries, if not handled properly, could have led to a global Ebola pandemic. There was also a separate Ebola outbreak affecting a remote region in the Democratic Republic of Congo (formerly Zaire), and it was only by November 21, 2014 that the outbreak was reported to have ended [2].

The Ebola Virus Disease (EVD), formally known as Ebola haemorrhagic fever and caused by the Ebola Virus, is very lethal with case fatalities ranging from 25% to 90%, with a mean of about 50% [2]. The 2014 EVD outbreak, though not the first but one of many other EVD outbreaks that have occurred in Africa since the first recorded outbreak of 1976, is the worst in terms of the numbers of Ebola cases and related deaths and the most complex [2]. About 9 months after the identification of a mysterious killer disease killing villagers in a small Guinean village as Ebola, the 2014 West African Ebola outbreak, as of December 24, 2014, had up to 19497 Ebola cases resulting in 7588 fatalities [1, 3, 4], a case fatality rate of about 38.9%. By December 2015, the number of Ebola Virus cases (including suspected, probable, and confirmed) stood at 28640 resulting in 11315 fatalities, a case fatality rate of 39.5% [3, 5].

Ebola Virus, the agent that causes EVD, is hypothesised to be introduced into the human population through contact with the blood, secretions, fluids from organs, and other body parts of dead or living animals infected with the virus (e.g., fruit bats, primates, and porcupines) [2, 6]. Human-to-human transmission can then occur through direct contact (via broken skin or mucous membranes such as eyes, nose, or mouth) with Ebola Virus infected blood, secretions, and fluids secreted through organs or other body parts, in, for example, saliva, vomit, urine, faeces, semen, sweat, and breast milk. Transmission can also be as a result of indirect contact with surfaces and materials, in, for example, bedding, clothing, and floor areas, or objects such as syringes, contaminated with the aforementioned fluids [2, 6].

When a healthy human (considered here to be susceptible) who has no Ebola Virus in them is exposed to the virus (directly or indirectly), the human may become infected, if transmission is successful. The risk of being infected with the Ebola Virus is (i) very low or not recognizable where there is casual contact with a feverish, ambulant, self-caring patient, for example, sharing the same public place, (ii) low where there is close face-to-face contact with a feverish and ambulant patient, for example, physical examination, (iii) high where there is close face-to-face contact without appropriate personal protective equipment (including eye protection) with a patient who is coughing or vomiting, has nosebleeds, or has diarrhea, and (iv) very high where there is percutaneous, needle stick, or mucosal exposure to virus-contaminated blood, body fluids, tissues, or laboratory specimens in severely ill or known positive patients. The symptoms of EVD may appear during the incubation period of the disease, estimated to be anywhere from 2 to 21 days [2, 7–9], with an estimated 8- to 10-day average incubation period, although studies estimate it at 9–11 days for the 2014 EVD outbreak in West Africa [10]. Studies have shown that, during the asymptomatic part of the Ebola Virus Disease, a human infected with the virus is not infectious and cannot transmit the virus. However, with the onset of symptoms, the human can transmit the virus and is hence infectious [2, 7]. The onset of symptoms commences the course of illness of the disease which can lead to death 6–16 days later [8, 9] or improvement of health with evidence of recovery 6–11 days later [8].

In the first few days of EVD illness (estimated at days 1–3 [11]), a symptomatic patient may exhibit symptoms common to those like the malaria disease or the flu (high fever, headache, muscle and joint pains, sore throat, and general weakness). Without effective disease management, between days 4 and 5 to 7, the patient progresses to gastrointestinal symptoms such as nausea, watery diarrhea, vomiting, and abdominal pain [10, 11]. Some or many of the other symptoms, such as low blood pressure, headaches, chest pain, shortness of breath, anemia, exhibition of a rash, abdominal pain, confusion, bleeding, and conjunctivitis, may develop [10, 11] in some patients. In the later phase of the course of the illness, days 7–10, patients may present with confusion and may exhibit signs of internal and/or visible bleeding, progressing towards coma, shock, and death [10, 11].

Recovery from EVD can be achieved, as evidenced by the less than 50% fatality rate for the 2014 EVD outbreak in West Africa. With no known cure, recovery is possible through effective disease management, the treatment of Ebola-related symptoms, and also the effective protection by the patient's immune response [7]. Some of the disease management strategies include hydrating patients by administering intravenous fluids and balancing electrolytes and maintaining the patient's blood pressure and oxygen levels. Other schemes used include blood transfusion (using an Ebola survivor's blood) and the use of experimental drugs on such patients (e.g., ZMAPP whose safety and efficacy have not yet been tested on humans). There are some other promising drugs/vaccines under trials [2]. Studies show that once a patient recovers from EVD they remain protected against the disease and are immune to it at least for a projected period because they develop antibodies that last for at least 10 years [7]. Once recovered, lifetime immunity is unknown or whether a recovered individual can be infected with another Ebola strain is unknown. However, after recovery, a person can potentially remain infectious as long as their blood and body fluids, including semen and breast milk, contain the virus. In particular, men can potentially transmit the virus through their seminal fluid, within the first 7 to 12 weeks after recovery from EVD [2]. Table 1 shows the estimated time frames and projected progression of the infection in an average EVD patient.

Table 1

A possible progression path of symptoms from  exposure to the Ebola Virus to treatment or death. Table shows a suggested transition and time frame in humans, of the virus, from exposure to incubation to symptoms development and recovery or death. This table is adapted based on the image in the Huffington post, via [11]. Superscript a: for the 2014 epidemic, the average incubation period is reported to be between 9 and 11 days [10]. Superscript b: other studies reported a mean of 4–10 days [8, 9].

Exposure	Incubation period	Course of illness			Recovery or death
	Range: 2 to 21 days from exposure	Range: 6 to 16 days from the end of the incubation period			Recovery: by the end of days 6–11	Death: by the end of days 6–16
		Probable	Early symptomatic	Late symptomatic
		Days 1–3	Days 4–7	Days 7–10
An individual comes in contact with an Ebola infected individual (dead or alive) or have been in the vicinity of someone who has been exposed.	Average of 8–11a days before symptoms are evident. Another estimate reports an average of 4–10b days.	Patients exhibit malaria-like or flu-like symptoms: for example, fever and weakness.	Patients progress to gastrointestinal symptoms: for example, nausea, watery diarrhea, vomiting, and abdominal pain. Other symptoms may include low blood pressure, anemia, headaches, chest pain, shortness of breath, exhibition of a rash, confusion, bleeding, and conjunctivitis.	Patients may present with confusion and may exhibit signs of internal and/or visible bleeding, potentially progressing towards coma, shock, and death.	Some patients may recover, while others will die.Recovery typically requires early intervention.

Given that there is no approved drug or vaccine out yet, local control of the Ebola Virus transmission requires a combined and coordinated control effort at the individual level, the community level, and the institutional/health/government level. Institutions and governments need to educate the public and raise awareness about risk factors, proper hand washing, proper handling of Ebola patients, quick reporting of suspected Ebola cases, safe burial practices, use of public transportation, and so forth. These education efforts need to be communicated with community/chief leaders who are trusted by members of the communities they serve. From a global perspective, a good surveillance and contact tracing program followed by isolation and monitoring of probable and suspected cases, with immediate commencement of disease management for patients exhibiting symptoms of EVD, is important if we must, in the future, elude a global epidemic and control of EVD transmission locally and globally [2]. It was by effective surveillance, contact tracing, and isolation and monitoring of probable and suspected cases followed by immediate supportive care for individuals and families exhibiting symptoms that the EVD was brought under control in Nigeria [17], Senegal, USA, and Spain [1].

Efficient control and management of any future EVD outbreaks can be achieved if new, more economical, and realizable methods are used to target and manage the dynamics of spread as well as the population sizes of those communities that may be exposed to any future Ebola Virus Disease outbreak. More realistic mathematical models can play a role in this regard, since analyses of such models can produce clear insight to vulnerable spots on the Ebola transmission chain where control efforts can be concentrated. Good models could also help in the identification of disease parameters that can possibly influence the size of the reproduction number of EVD. Existing mathematical models for Ebola [14, 16, 18–21] have been very instrumental in providing mathematical insight into the dynamics of Ebola Virus transmission. Many of these models have also been helpful in that they have provided methods to derive estimates for the reproduction number for Ebola based on data from the previous outbreaks. However, few of the models have taken into account the fact that institution of quarantine states or treatment centres will affect the course of the epidemic in the population [16]. It is our understanding that the way the disease will spread will be determined by the initial and continual response of the health services in the event of the discovery of an Ebola disease case. The objective of this paper is to derive a comprehensive mathematical model for the dynamics of Ebola transmission taking into consideration what is currently known of the disease. The primary objective is to derive a formula for the reproduction number for Ebola Virus Disease transmission in view of providing a more complete and measurable index for the spread of Ebola and to investigate the level of impact of surveillance, contact tracing, isolation, and monitoring of suspected cases, in curbing disease transmission. The model is formulated in a way that it is extendable, with appropriate modifications, to other disease outbreaks with similar characteristics to Ebola, requiring such contact tracing strategies. Our model differs from other mathematical models that have been used to study the Ebola disease [14, 15, 18, 20–22] in that it captures the quarantined Ebola Virus Disease patients and provides possibilities for those who escaped quarantine at the onset of the disease to enter quarantine at later stages. To the best of our knowledge, this is the first integrated ordinary differential equation model for this kind of communicable disease of humans. Our final result would be a formula for the basic reproduction number of Ebola that depends on the disease parameters.

The rest of the paper is divided up as follows. In Section 2, we outline the derivation of the model showing the state variables and parameters used and how they relate together in a conceptual framework. In Section 3, we present a mathematical analysis of the derived model to ascertain that the results are physically realizable. We then reparameterise the model and investigate the existence and linear stability of steady state solution, calculate the basic reproduction number, and present some special cases. In Section 4, we present a discussion on the parameters of the model. In Section 5, we carry out some numerical simulations based on the selected feasible parameters for the system and then round up the paper with a discussion and conclusion in Section 6.

2. The Mathematical Model

2.1. Description of Model Variables

We divide the human population into 11 states representing disease status and quarantine state. At any time t there are the following.

(1) Susceptible Individuals. Denoted by S, this class also includes false probable cases, that is, all those individuals who would have displayed early Ebola-like symptoms but who eventually return a negative test for Ebola Virus infection.

(2) Suspected Ebola Cases. The class of suspected EVD patients comprises those who have come in contact with, or been in the vicinity of, anybody who is known to have been sick or died of Ebola. Individuals in this class may or may not show symptoms. Two types of suspected cases are included: the quarantined suspected cases, denoted by S , and the nonquarantined suspected case, denoted by S . Thus a suspected case is either quarantined or not.

(3) Probable Cases. The class of probable cases comprises all those persons who at some point were considered suspected cases and who now present with fever and at least three other early Ebola-like symptoms. Two types of probable cases are included: the quarantined probable cases, denoted by P , and the nonquarantined probable cases, P . Thus a probable case is either quarantined or not. Since the early Ebola-like symptoms of high fever, headache, muscle and joint pains, sore throat, and general weakness can also be a result of other infectious diseases such as malaria or flu, we cannot be certain at this stage whether or not the persons concerned have Ebola infection. However, since the class of probable persons is derived from suspected cases, and to remove the uncertainties, we will assume that probable cases may eventually turn out to be EVD patients and if that were to be the case, since they are already exhibiting some symptoms, they can be assumed to be mildly infectious.

(4) Confirmed Early Symptomatic Cases. The class of confirmed early asymptomatic cases comprises all those persons who at some point were considered probable cases and a confirmatory laboratory test has been conducted to confirm that there is indeed an infection with Ebola Virus. This class is called confirmed early symptomatic because all that they have as symptoms are the early Ebola-like symptoms of high fever, headache, muscle and joint pains, sore throat, and general weakness. Two types of confirmed early symptomatic cases are included: the quarantined confirmed early symptomatic cases C and the nonquarantined confirmed early symptomatic cases C . Thus a confirmed early symptomatic case is either quarantined or not. The class of confirmed early symptomatic individuals may not be very infectious.

(5) Confirmed Late Symptomatic Cases. The class of confirmed late symptomatic cases comprises all those persons who at some point were considered confirmed early symptomatic cases and in addition the persons who now present with most or all of the later Ebola-like symptoms of vomiting, diarrhea, stomach pain, skin rash, red eyes, hiccups, internal bleeding, and external bleeding. Two types of confirmed late symptomatic cases are included: the quarantined confirmed late symptomatic cases I and the nonquarantined confirmed late symptomatic cases I . Thus a confirmed late symptomatic case is either quarantined or not. The class of confirmed late symptomatic individuals may be very infectious and any bodily secretions from this class of persons can be infectious to other humans.

(6) Removed Individuals. Three types of removals are considered, but only two are related to EVD. The removals related to the EVD are confirmed individuals removed from the system through disease induced death, denoted by R , or confirmed cases that recover from the infection denoted by R . Now, it is known that unburied bodies or not yet cremated cadavers of EVD victims can infect other susceptible humans upon contact [7]. Therefore, the cycle of infection really stops only when a cadaver is properly buried or cremated. Thus members from class, R , representing dead bodies or cadavers of EVD victims are considered removed from the infection chain, and consequently from the system, only when they have been properly disposed of. The class, R , of individuals who beat the odds and recover from their infection are considered removed because recovery is accompanied with the acquisition of immunity so that this class of individuals are then protected against further infection [7] and they no longer join the class of susceptible individuals. The third type of removal is obtained by considering individuals who die naturally or due to other causes other than EVD. These individuals are counted as R .

The state variables are summarized in Notations.

2.2. The Mathematical Model

A compartmental framework is used to model the possible spread of EVD within a population. The model accounts for contact tracing and quarantining, in which individuals who have come in contact or have been associated with Ebola infected or Ebola-deceased humans are sought and quarantined. They are monitored for twenty-one days during which they may exhibit signs and symptoms of the Ebola Virus or are cleared and declared free. We assume that most of the quarantining occurs at designated makeshift, temporal, or permanent health facilities. However, it has been documented that others do not get quarantined, because of fear of dying without a loved one near them or fear that if quarantined they may instead get infected at the centre, as well as traditional practices and belief systems [14, 16, 22]. Thus, there may be many within communities who remain nonquarantined, and we consider these groups in our model. In all the living classes discussed, we will assume that natural death, or death due to other causes, occurs at constant rate μ where 1/μ is approximately the life span of the human.

2.2.1. The Susceptible Individuals

The number of susceptible individuals in the population decreases when this population is exposed by having come in contact with or being associated with any of the possibly infectious cases, namely, infected probable case, confirmed case, or the cadaver of a confirmed case. The density increases when some false suspected individuals (a proportion of 1 − θ 2 of nonquarantined and 1 − θ 6 of quarantined) and probable cases (a proportion of 1 − θ 3 of nonquarantined individuals and 1 − θ 7 of quarantined individuals) are eliminated from the suspected and probable case list. We also assume a constant recruitment rate Π as well as natural death, or death due to other causes. Therefore the equation governing the rate of change with time within the class of susceptible individuals may be written aswhere λ is the force of infection and the rest of the parameters are positive and are defined in Notations. We identify two types of total populations at any time t: (i) the total living population, H , and (ii) the total living population including the cadavers of Ebola Virus victims that can take part in the spread of EVD, H. Thus at each time t we have Since the cadavers of EVD victims that have not been properly disposed of are very infectious, the force of infection must then also take this fact into consideration and be weighted with H instead of H . The force of infection takes the following form:where H > 0 is defined above and the parameters ρ , ρ , τ , τ , ξ , ξ , and a are positive constants as defined in Notations. There are no contributions to the force of infection from the R class because it is assumed that once a person recovers from EVD infection, the recovered individual acquires immunity to subsequent infection with the same strain of the virus. Although studies have suggested that recovered men can potentially transmit the Ebola Virus through seminal fluids within the first 7–12 weeks of recovery [2], and mothers through breast milk, we assume, here, that, with education, survivors who recover would have enough information to practice safe sexual and/or feeding habits to protect their loved ones until completely clear. Thus recovered individuals are considered not to contribute to the force of infection.

2.2.2. The Suspected Individuals

A fraction 1 − θ 1 of the exposed susceptible individuals get quarantined while the remaining fraction are not. Also, a fraction θ 2 (resp., θ 6) of the nonquarantined (resp., quarantined) suspected individuals become probable cases at rate α while the remainder 1 − θ 2 (resp., 1 − θ 6) do not develop into probable cases and return to the susceptible pool. For the quarantined individuals, we assume that they are being monitored, while the suspected nonquarantined individuals are not. However, as they progress to probable cases (at rates α and α ), a fraction θ 2 of these humans will seek the health care services as symptoms commence and become quarantined while the remainder θ 2 remain nonquarantined. Thus the equation governing the rate of change within the two classes of suspected persons then takes the following form:

In the context of this model we make the assumption that once quarantined, the individuals stay quarantined until clearance and are released, or they die of the infection. Notice that θ 2 = θ 2 + θ 2, so that 1 − θ 2 + θ 2 + θ 2 = 1.

2.2.3. The Probable Cases

The fractions θ 2 and θ 6 of suspected cases that become probable cases increase the number of individuals in the probable case class. The population of probable cases is reduced (at rates β and β ) when some of these are confirmed to have the Ebola Virus through laboratory tests at rates α and α . For some, proportions 1 − θ 3 and 1 − θ 7, the laboratory tests are negative and the probable individuals revert to the susceptible class. From the proportion θ 3 of nonquarantined probable cases whose tests are positive for the Ebola Virus (i.e., confirmed for EVD), a fraction, θ 3, become quarantined while the remainder, θ 3, remain nonquarantined. So θ 3 = θ 3 + θ 3. Thus the equation governing the rate of change within the classes of probable cases takes the following form:

2.2.4. The Confirmed Early Symptomatic Cases

The fractions θ 3 and θ 7 of the probable cases become confirmed early symptomatic cases thus increasing the number of confirmed cases with early symptoms. The population of early symptomatic individuals is reduced when some recover at rates r for the nonquarantined cases and r for the quarantined cases. Others may see their condition worsening and progress and become late symptomatic individuals, in which case they enter the full blown late symptomatic stages of the disease. We assume that this progression occurs at rates γ or γ , respectively, which are the reciprocal of the mean time it takes for the immune system to either be completely overwhelmed by the virus or be kept in check via supportive mechanism. A fraction 1 − θ 4 of the confirmed nonquarantined early symptomatic cases will be quarantined as they become late symptomatic cases, while the remaining fraction θ 4 escape quarantine due to lack of hospital space or fear and belief customs [16, 22] but become confirmed late symptomatic cases in the community. Thus the equation governing the rate of change within the two classes of confirmed early symptomatic cases takes the following form:

2.2.5. The Confirmed Late Symptomatic Cases

The fractions θ 4 and θ 8 of confirmed early symptomatic cases who progress to the late symptomatic stage increase the number of confirmed late symptomatic cases. The population of late symptomatic individuals is reduced when some of these individuals are removed. Removal could be as a result of recovery at rates proportional to r and r or as a result of death because the EVD patient's conditions worsen and the Ebola Virus kills them. The death rates are assumed proportional to δ and δ . Additionally, as a control effort or a desperate means towards survival, some of the nonquarantined late symptomatic cases are removed and quarantined at rate σ . In our model, we assume that Ebola-related death only occurs at the late symptomatic stage. Additionally, we assume that the confirmed late symptomatic individuals who are eventually put into quarantine at this late period (removing them from the community) may not have long to live but may have a slightly higher chance at recovery than when in the community and nonquarantined. Since recovery confers immunity against the particular strain of the Ebola Virus, individuals who recover become refractory to further infection and hence are removed from the population of susceptible individuals. Thus the equation governing the rate of change within the two classes of confirmed late symptomatic cases takes the following form:

2.2.6. The Cadavers and the Recovered Persons

The dead bodies of EVD victims are still very infectious and can still infect susceptible individuals upon effective contact [2]. Disease induced deaths from the class of confirmed late symptomatic individuals occur at rates δ and δ and the cadavers are disposed of via burial or cremation at rate b. The recovered class contains all individuals who recover from EVD. Since recovery is assumed to confer immunity against the 2014 strain (the Zaire Virus) [7] of the Ebola Virus, once an individual recovers, they become removed from the population of susceptible individuals. Thus the equation governing the rate of change within the two classes of recovered persons and cadavers takes the following form:

The population of humans who die either naturally or due to other causes is represented by the variable R and keeps track of all natural deaths, occurring at rate μ, from all the living population classes. This is a collection class. Another collection class is the class of disposed Ebola-related cadavers, disposed at rate b. These collection classes satisfy the equations

Putting all the equations together we have where and all other parameters and state variables are as in Notations.

Suitable initial conditions are needed to completely specify the problem under consideration. We can, for example, assume that we have a completely susceptible population, and a number of infectious persons are introduced into the population at some point. We can, for example, have that

Class D is used to keep count of all the dead that are properly disposed of, class R is used to keep count of all the deaths due to EVD, and class R is used to keep count of the deaths due to causes other than EVD infection. The rate of change equation for the two groups of total populations is obtained by using (2) and (3) and adding up the relevant equations from (11) to (23) to obtain where H is the total living population and H is the augmented total population adjusted to account for nondisposed cadavers that are known to be very infectious. On the other hand, if we keep count of all classes by adding up (11)–(23), the total human population (living and dead) will be constant if Π = 0. In what follows, we will use the classes R , R , and D , comprising classes of already dead persons, only as place holders, and study the problem containing the living humans and their possible interactions with cadavers of EVD victims as often is the case in some cultures in Africa, and so we cannot have a constant total population. Note that (26) can also be written as dH/dt = Π − μH − bR .

2.2.7. Infectivity of Persons Infected with EVD

Ebola is a highly infectious disease and person to person transmission is possible whenever a susceptible person comes in contact with bodily fluids from an individual infected with the Ebola Virus. We therefore define effective contact here generally to mean contact with these fluids. The level of infectivity of an infected person usually increases with duration of the infection and severity of symptoms and the cadavers of EVD victims are the most infectious [23]. Thus we will assume in this paper that probable persons who indeed are infected with the Ebola Virus are the least infectious while confirmed late symptomatic cases are very infectious and the level of infectivity will culminate with that of the cadaver of an EVD victim. While under quarantine, it is assumed that contact between the persons in quarantine and the susceptible individuals is minimal. Thus though the potential infectivity of the corresponding class of persons in quarantine increases with disease progression, their effective transmission to members of the public is small compared to that from the nonquarantined class. It is therefore reasonable to assume that any transmission from persons under quarantine will affect mostly health care providers and use that branch of the dynamics to study the effect of the transmission of the infections to health care providers who are here considered part of the total population. In what follows we do not explicitly single out the infectivity of those in quarantine but study general dynamics as derived by the current modelling exercise.

3. Mathematical Analysis

3.1. Well-Posedness, Positivity, and Boundedness of Solution

In this subsection we discuss important properties of the model such as well-posedness, positivity, and boundedness of the solutions. We start by defining what we mean by a realistic solution.

Definition 1 (realistic solution).

A solution of system (27) or equivalently system comprising (11)–(21) is called realistic if it is nonnegative and bounded.

It is evident that a solution satisfying Definition 1 is physically realizable in the sense that its values can be measured through data collection. For notational simplicity, we use vector notation as follows: let x = (S, S , S , P , P , C , C , I , I , R , R ) be a column vector in ℝ 11 containing the 11 state variables, so that, in this notation, x 1 = S, x 2 = S ,…, x 11 = R . Let f(x) = (f 1(x), f 2(x),…, f 11(x)) be the vector valued function defined in ℝ 11 so that in this notation f 1(x) is the right-hand side of the differential equation for first variable S, f 2(x) is the right-hand side of the equation for the second variable x 2 = S ,…, and f 11(x) is the right side of the differential equation for the 11th variable x 11 = R , and so is precisely system (11)–(21) in that order with prototype initial conditions (24). We then write the system in the formwhere x : [0, ∞) → ℝ 11 is a column vector of state variables and f : ℝ 11 → ℝ 11 is the vector containing the right-hand sides of each of the state variables as derived from corresponding equations in (11)–(21). We can then have the following result.

Lemma 2 .

The function f in (27) is Lipschitz continuous in x.

Proof

Since all the terms in the right-hand side are linear polynomials or rational functions of nonvanishing polynomial functions, and since the state variables, S, S , P , C , I , and R , are continuously differentiable functions of t, the components of the vector valued function f of (27) are all continuously differentiable. Further, let L(x, y; θ) = {x + θ(y − x) : 0 ≤ θ ≤ 1}. Then L(x, y; θ) is a line segment that joins points x to the point y as θ ranges on the interval [0,1]. We apply the mean value theorem to see that where f′(z; y − x) is the directional derivative of the function f at the mean value point z in the direction of the vector y − x. But, where e is the kth coordinate unit vector in ℝ + 11. It is now a straightforward computation to verify that since ℝ + 11 is a convex set, and taking into consideration the nature of the functions f , i = 1,…, 11, all the partial derivatives are bounded and so there exist M > 0 such that and so there exist M > 0 such that and hence f is Lipschitz continuous.

Theorem 3 (uniqueness of solutions).

The differential equation (27) has a unique solution.

By Lemma 2, the right-hand side of (27) is Lipschitzian; hence a unique solution exists by existence and uniqueness theorem of Picard. See, for example, [24].

Theorem 4 (positivity).

The region ℝ + 11 wherein solutions defined by (11)–(21) are defined is positively invariant under the flow defined by that system.

We show that each trajectory of the system starting in ℝ + 11 will remain in ℝ + 11. Assume for a contradiction that there exists a point t 1 ∈ [0, ∞) such that S(t 1) = 0, S′(t 1) < 0 (where the prime denotes differentiation with respect to time) but for 0 < t < t 1, S(t) > 0, and S (t) > 0, S (t) > 0, P (t) > 0, P (t) > 0, C (t) > 0, C (t) > 0, I (t) > 0, I (t) > 0, R (t) > 0, and R (t) > 0. So, at the point t = t 1, S(t) is decreasing from the value zero in which case it will go negative. If such an S will satisfy the given differential equation, then we haveand so S′(t 1) > 0 contradicting the assumption that S′(t 1) < 0. So no such t 1 exist. The same argument can be made for all the state variables. It is now a simple matter to verify, using techniques as explained in [24], that whenever we start system (27), with nonnegative initial data in ℝ + 11, the solution will remain nonnegative for all t > 0 and that if x 0 = 0, the solution will remain x = 0  ∀t > 0, and the region ℝ + 11 is indeed positively invariant.

The last two theorems have established the fact that, from a mathematical and physical standpoint, the differential equation (27) is well-posed. We next show that the nonnegative unique solutions postulated by Theorem 3 are indeed realistic in the sense of Definition 1.

Theorem 5 (boundedness).

The nonnegative solutions characterized by Theorems 3 and 4 are bounded.

It suffices to prove that the total living population size is bounded for all t > 0. We show that the solutions lie in the bounded regionFrom the definition of H given in (2), if H is bounded, the rest of the state variables that add up to H will also be bounded. From (25) we haveThus, from (34), we see that, whatever the size of H (0), H (t) is bounded above by a quantity that converges to Π/μ as t → ∞. In particular, if μH (0) < Π, then H (t) is bounded above by Π/μ, and for all initial conditions Thus H (t) is nonnegative and bounded.

Remark 6 .

Starting from the premise that H (t) ≥ 0 for all t > 0, Theorem 5 establishes boundedness for the total living population and thus by extension verifies the positive invariance of the positive octant in ℝ 11 as postulated by Theorem 4, since each of the variables functions S, S , P , C , I , and R , where ∗ ∈ {N, Q, R}, is a subset of H .

3.2. Reparameterisation and Nondimensionalisation

The only physical dimension in our system is that of time. But we have state variables which depend on the density of humans and parameters which depend on the interactions between the different classes of humans. A state variable or parameter that measures the number of individuals of certain type has dimension-like quantity associated with it [25]. To remove the dimension-like character on the parameters and variables, we make the following change of variables:where We then define the dimensionless parameter groupingsThe force of infection λ then takes the formThis leads to the equivalent system of equations and the total populations satisfy the scaled equation where b 8 > 0 if it is assumed that the rate of disposal of Ebola Virus Disease victims, b, is larger than the natural human death rate, μ. The scaled or dimensionless parameters are then as follows:

3.3. The Steady State Solutions and Linear Stability

The steady state of the system is obtained by setting the right-hand side of the scaled system to zero and solving for the scalar equations. Let x = (s , s , s , p , p , c , c , i , i , r , r , h ) be a steady state solution of the system. Then, (43), (45), and (47) indicate thatand we can use any of these as a parameter to derive the values of the other steady state variables. We use the variables p and s as parameters to obtain the expressionswhereHere the solution for the scaled total living (h ) and scaled living and Ebola-deceased (h) populations is, respectively, obtained by equating the right-hand sides of (53) and (54) to zero, while that for s is obtained by adding up (40), (41), and (42). It is easy to verify from reparameterisation (38) that the parameter groupings B 3 and B 4 are both nonnegative. In fact, showing that B 3 > 0 and B 4 > 0.

To obtain a value for p and s , we substitute all computed steady state values, (56) and (57), into (41) and (42). The expression for λ s in terms of p and s is obtained from (39). Performing the aforementioned procedures leads to the two equations whereNext, we solve (60) and (61) simultaneously, which clearly differ in some of their coefficients, to obtain the expressions for p and s . Quickly observe that the two equations may be reduced to one such thatTwo possibilities arise: either (i) 1 − B 1 p − B 2 s = 0 or (ii) α p (1 − θ 1) − α s θ 1 = 0. The first condition leads to the systemHowever, the two equations are equivalent to h = 0 and s = 0 (see (57)), which are unrealistic, based on our constant population recruitment model. Hence, we only consider the second possibility, which yields the relationso that substituting (65) into (60) yieldswhere with

Remark 7 .

It can be shown that B 2 < B 4. In fact,

Thus, if (1 − θ 1)α (B 4 − B 2) > θ 1 α (B 1 − B 3), then z > 1. This will hold if B 3 > B 1. In the case where B 3 < B 1, we will require that B 4 − B 2 be greater than (θ 1 α /(1 − θ 1)α )(B 1 − B 3).

We identify R 0 as the unique threshold parameter of the system as follows.

Lemma 8 .

The parameter R 0 defined in (69) is the unique threshold parameter of the system whenever z > 1.

If z > 1, then ℛ = zR 0 > 1 whenever R 0 > 1 and the existence or nonexistence of a realistic solution of the form of (66) is determined solely by the size of R 0.

The rest of the steady states are then obtained by using these values for p and s given by (66) in (65) and (57) to obtain the following:where x, y, and z are as defined in (68). We have proved the following result.

Theorem 9 (on the existence of equilibrium solutions).

System (40)–(52) has at least two equilibrium solutions: the disease-free equilibrium x = E = (1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1) and an endemic equilibrium x = E = (s , s , s , p , p , c , c , i , i , r , r , h ). The endemic equilibrium, E , exists and is realistic only when the threshold parameters R 0 and ℛ, given by (69), are of appropriate magnitude.

The stability of the steady states is governed by the sign of the eigenvalues of the linearizing matrix near the steady state solutions. If J(x ) is the Jacobian matrix at the steady state x , then we havewhere . Thus if ζ is an eigenvalue of the linearized system at the disease-free state, then ζ is obtained by the solvability conditionan equation involving a polynomial of degree 12 in ζ, where P 9(ζ) is a polynomial of degree 9 in ζ, given byNow, all we need to know at this stage is whether there is solution of (73) for ζ with positive real part which will then indicate the existence of unstable perturbations in the linear regime. The coefficients of polynomial (73) can give us vital information about the stability or instability of the disease-free equilibrium. For example, by Descartes' rule of signs, a sign change in the sequence of coefficients indicates the presence of a positive real root which in the linear regime signifies the presence of exponentially growing perturbations. We can write polynomial equation (73) in the formwhereand we can see that c 0 changes sign from positive to negative when R 0 increases from values of R 0 < 1 through R 0 = 1 to values of R 0 > 1 indicating a change in stability of the disease-free equilibrium as R 0 increases from unity.

3.4. The Basic Reproduction Number

A threshold parameter that is of essential importance to infectious disease transmission is the basic reproduction number denoted by R 0. R 0 measures the average number of secondary clinical cases of infection generated in an absolutely susceptible population by a single infectious individual throughout the period within which the individual is infectious [26-29]. Generally, the disease eventually disappears from the community if R 0 < 1 (and in some situations there is the occurrence of backward bifurcation) and may possibly establish itself within the community if R 0 > 1. The critical case R 0 = 1 represents the situation in which the disease reproduces itself thereby leaving the community with a similar number of infection cases at any time. The definition of R 0 specifically requires that initially everybody but the infectious individual in the population be susceptible. Thus, this definition breaks down within a population in which some of the individuals are already infected or immune to the disease under consideration. In such a case, the notion of reproduction number ℜ becomes useful. Unlike R 0 which is fixed, ℜ may vary considerably with disease progression. However, ℜ is bounded from above by R 0 and it is computed at different points depending on the number of infected or immune cases in the population.

One way of calculating R 0 is to determine a threshold condition for which endemic steady state solutions to the system under study exist (as we did to derive (69)) or for which the disease-free steady state is unstable. Another method is the next-generation approach where R 0 is the spectral radius of the next-generation matrix [26]. Using the next-generation approach, we identify all state variables for the infection process, p , p , c , c , i , i , r , r , and h. The transitions from s , s to p , p are not considered new infections but rather a progression of the infected individuals through the different stages of disease compartments. Hence, we identify terms representing new infections from the above equations and rewrite the system as the difference of two vectors and , where consists of all new infections and consists of the remaining terms or transitions between states. That is, we set , where x is the vector of state variables corresponding to new infections: x = (s, s , s , p , p , c , c , i , i , r , r , h). This gives rise towhere the force of infection λ is given by (39). To obtain the next-generation operator, FV −1, we must calculate and evaluated at the disease-free equilibrium position, where s = 1 = h, s = s = p = p = c = c = i = i = r = r = 0. The basic reproduction number is then the spectral radius of the next-generation matrix FV −1. Thus if ϱ(FV −1) is the spectral radius of the matrix FV −1, thenas computed before.

The expression for R 0 has two parts. The first part measures the number of new EVD cases generated by an infected nonquarantined human. It is the product of θ 1 (the proportion of susceptible individuals who become suspected but remain nonquarantined), B 5 (which indicates the contacts from this proportion of individuals with infected individuals at various stages of the disease), and 1/α (which is the average duration a human remains as a suspected nonquarantined individual). In the same way, the second part can be interpreted likewise.

The stability of the endemic steady state is obtained by calculating the eigenvalues of the linearized matrix evaluated at the endemic state. The computations soon become very complicated because of the size of the system and we proceed with a simplification of the system.

3.5. Pseudo-Steady State Approximation

Ebola Virus Disease is a very deadly infection that normally kills most of its victims within about 21 days of exposure to the infection. Thus when compared with the life span of the human, the elapsed time representing the progression of the infection from first exposure to death is short when compared to the total time required as the life span of the human. Thus we set μ ≈ 1/life  span  of  human, so that the rates α , β , and so forth will all be such that 1/rate ≈ resident  time  in  given  state, some of which will be short compared with the life span of the human. It is therefore reasonable to assume thatso that the scaling above renders some of the state variables essentially at equilibrium. That is, the quantities, 1/β , 1/β , 1/γ , 1/γ , 1/δ , 1/δ , and b/μ, may be regarded as small parameters so that, in the corresponding equations (43)–(48) and (50), the state variables modelled by these equations are essentially in equilibrium and we can evoke the Michaelis-Menten pseudo-steady state hypothesis [30]. To proceed, we make the pseudoequilibrium approximationto have the reduced system and the total population and the force of infection also reduce accordingly. In particular, we havewhere the variables s and the augmented population h satisfy the differential equations (81) and (85), respectively.

System (81)–(85) has the same steady states solutions as the original system if we combine it with (80). However, on its own, it represents a pseudo-steady state approximation [30] of the original system. Clearly the reduced system has two realistic steady states: E dfe and E , so that if E = (s , s , s , r , h ) is a steady state solution, thenwhere, following the same method as was done in the full system,All coefficients are as defined in (58). When steady states (88) are rendered in parameters of the reduced system, taking into consideration the fact that, from (68), z = B 3 y + B 4 x and B 1 y + B 2 x = 1, we get

The stability of the steady states is determined by the eigenvalues of the linearized matrix of the reduced system evaluated at the steady state x = (s , s , s , r , h ). If J(x ) is the Jacobian of the system near the steady state x , thenwhere and The asterisk is used to indicate that the quantities so calculated are evaluated at the steady state. We can perform a stability analysis on the reduced system by noting that if ζ is an eigenvalue of (90), then ζ satisfies the polynomial equationwhere Now, the signs of the zeros of (92) will depend on the signs of the coefficients Q , i ∈ {0,1, 2}. We now examine these.

At the disease-free state where s = 1 = h , s = s = r = 0, or equivalently R 0 = 1, we have x = x dfe = (1,0, 0,0, 1) so that C 3(x dfe) = 0, C 4(x dfe) = B 5, C 5(x dfe) = B 6, C 6(x dfe) = 0, and (92) becomeswhere The roots of (94) are −1, −1, −1, and , showing that there is one positive real solution as R 0 increases beyond unity and the disease-free equilibrium loses stability at R 0 = 1. For the local stability when R 0 ≤ 1, the additional requirement Q dfe > 0 is necessary.

At the endemic steady state, and in the original scaled parameter groupings of the system, x = x = (s , s , s , r , h ), the coefficients of (92) simplify accordingly and we havewherewhere Now the necessary and sufficient conditions that will guarantee the stability of the nontrivial steady state x will be the Routh-Hurwitz criteria which, in the present parameterizations, areWith this characterization, we can then explore special cases of intervention.

3.6. Some Special Cases

All initial suspected cases are quarantined: that is, θ 1 = 0. In this case we see that s = 0 and we have only the right branch of our flow chart in Figure 1. We have here a problem involving infections only at the treatment centres. Mathematically, we then haveand (96) becomesshowing that all solutions of the equation P 5(ζ, x ) = 0 are negative or have negative real parts whenever they are complex, indicating that the nontrivial steady state is stable to small perturbations whenever R 0 > 1. In this case we can regard an increase in R 0 as an increase in the parameter grouping B 6.

Figure 1

Conceptual framework showing the relationships between the different compartments that make up the different population of individuals and actors in the case of an EVD outbreak. Susceptible individuals include false suspected and probable cases. True suspects and probable cases are confirmed by a laboratory test and the confirmed cases can later develop symptoms and die of the infection or recover to become immune to the infection. Humans can also die naturally or due to other causes. Nonquarantined cases can become quarantined through intervention strategies. Others run the course of the illness from infection to death without being quarantined. Flow from compartment to compartment is as explained in the text.

All initial suspected cases escape quarantine: that is, θ 1 = 1. In this case we see that s = 0 and initially we will be on the left branch of our flow chart in Figure 1. The strength of the present model is that, based on its derivation, it is possible for some individuals to eventually enter quarantine as the systems wake up from sleep and control measures kick into place. Mathematically, we then haveand (96) becomes again showing that all solutions of the equation P 5(ζ, x ) = 0 are negative or have negative real parts whenever they are complex, indicating that the nontrivial steady state is stable to small perturbations whenever R 0 > 1. In this case we can regard an increase in R 0 as an increase in the parameter grouping B 5. R 0 in this case appears larger than in the previous cases.

The rate at which suspected individuals become probable cases is the same: that is, α = α . In this case the flow from being a suspected case to a probable case is the same in all circumstances, irrespective of whether or not one is quarantined or nonquarantined. Mathematically, we have that α = α and (96) becomesIn this case as well, all solutions of the equation P 5(ζ, x ) = 0 either are negative or have negative real parts whenever R 0 > 1 and z > 1, showing that in this case again the steady state is stable to small perturbations. In this particular case, an increase in R 0 can be regarded as an increase in the two parameter groupings B 5 and B 6.

4. Parameter Discussion

Some parameter values were chosen based on estimates in [15, 20], on the 2014 Ebola outbreak, while others were selected from past estimates (see [9, 14, 18]) and are summarized in Table 2. In [20], an estimate based on data primarily from March to August 20th yielded the following average transmission rates and 95% confidence intervals: 0.27 (0.27,0.27) per day in Guinea, 0.45 (0.43,0.48) per day in Sierra Leone, and 0.28 (0.28,0.29) per day in Liberia. In [15], the number of cases of the 2014 Ebola outbreak data (up until early October) was fitted to a discrete mathematical model, yielding estimates for the contact rates (per day) in the community and hospital (considered quarantined) settings as 0.128 in Sierra Leone for nonquarantined cases (and 0.080 for quarantined cases, about a 61% reduction) while the rates in Liberia were 0.160 for nonquarantined cases (and 0.062 for quarantined cases, close to a 37.5% drop). The models in [15, 20] did not separate transmission based on early or late symptomatic EVD cases, which was considered in our model. Based on the information, we will assume that effective contacts between susceptible humans and late symptomatic EVD patients in the communities fall in the range [0.12,0.48] per day, which contains the range cited in [16]. Thus ξ ∈ [0.12,0.48]. However, we will consider scenarios in which this parameter varies. Furthermore, if we assume about a 37.5% to 61% reduction in effective contact rates in the quarantined settings, then we can assume that ξ = ϕ ξ , where ϕ ∈ [0.0375,0.062]. However, to understand how effective quarantining Ebola patients is, general values of ϕ ∈ [0,1] can be considered. Under these assumptions, a small value of ϕ will indicate that quarantining was effective, while values of ϕ close to 1 will indicate that quarantining patients had no effect in minimizing contacts and reducing transmissions.

Table 2

Parameters, baseline values, and ranges of baseline values with references.

Parameters	Baseline values	Range of values	Reference
Π	3, 555	Varies
ρ N	Varies	Varies
ρ Q	Varies	Varies
ξ N	0.27	[0.12,0.48]	Estimated∗
τ N	q N ξ N, q N = 0.75	q N ∈ [0,1]	Variable
ξ Q	ϕ Q ξ N, ϕ Q = 0.5	ϕ Q ∈ [0,1]	Estimated∗
τ Q	q Q ξ Q, q Q = 0.75	q Q ∈ [0,1]	Variable
θ i	θ i ∈ [0,1], i = 1,2,…, 7	θ i ∈ [0,1], i = 1,2,…, 7	Variable
α N	1/10	[1/11,1/4]	[8–10]
α Q	1/10	[1/11,1/4]	[8–10]
β N	0.5	[1/3,1]	[10]
β Q	1/2	[1/3,1]	[10]
γ N	1/3	[1/4,1/2]	[10, 12]
γ Q	1/3	[1/4,1/2]	[10, 12]
δ N	1/3	[1/4,1/2]	[10, 12]
δ Q	1/3	[1/4,1/2]	[10, 12]
μ	1/(60 × 365)	[1/(40 × 365), 1/(70 × 365)] day−1	[13]
b	1/2.5	[1/4.50,1/2] day−1	[14, 15]
a D	0.3000	[0.111,0.489] day−1	[15]
r EQ	0.5366	[0.4829,0.5903]	[16]
r EN	ωr EQ, ω = 0.88	ω ∈ [0,1]	Estimate
r LN	κr EN, κ = 0.02	κ ≪ 1	Estimate
r LQ	κr EQ, κ = 0.02	κ ≪ 1	Estimate
σ N	0.5	[1/3,1]	Estimate

Estimates discussed in Section 4.

Since patients with EVD at the onset of symptoms are less infectious than EVD patients in the later stages of symptoms [2, 10], we assume that the effective contact rate between confirmed nonquarantined early symptomatic individuals and susceptible individuals, denoted by τ , is proportional to ξ with proportionality constant q ∈ [0,1]. Likewise, we assume that the effective contact rate between confirmed quarantined early symptomatic individuals and susceptible individuals is proportional to ξ with proportionality constant q ∈ [0,1]. Thus, τ = q ξ and τ = q ξ , with 0 < q , q < 1.

The range, for the parameter a , of the effective contact rate between cadavers of confirmed late symptomatic individuals and susceptible individuals was chosen to be [0.111,0.489] per day, where 0.111 is the rate estimated in [15] for Sierra Leone and 0.489 is that for Liberia. With control measures and education in place, these rates can be much lower.

The incubation period of EVD is estimated to be between 2 and 21 days [2, 7–9], with a mean of 4–10 days reported in [8, 9]. In [10], a mean incubation period of 9–11 days was reported for the 2014 EVD. Here, we will consider a range from 4 to 11 days, with a mean of about 10 days used as the baseline value. Thus we will consider that α and α are in the range [1/11,1/4]. At the end of the incubation period, early symptoms may emerge 1–3 days later [10], with a mean of 2 days. Thus the mean rate at which nonquarantined (β ) and quarantined (β ) suspected cases become probable cases lies in [1/3,1] [12]. About 1 or 2 to 4 days later after the early symptoms, more severe symptoms may develop so that the rates at which nonquarantined (γ ) and quarantined (γ ) probable cases become confirmed cases lie in [1/4,1/2].

The parameter γ measures the rate at which early symptomatic individuals leave that class. This could be as a result of recovery or due to increase and spread of the virus within the human. It takes about 2 to 4 days to progress from the early symptomatic stage to the late symptomatic stage, so that γ , γ ∈ [1/4,1/2], which can be assumed to be the reciprocal of the mean time it takes from when the immune system is either completely overwhelmed by the virus or kept in check via supportive mechanism. Severe symptoms are followed either by death after about an average of two to four days beyond entering the late symptomatic stage or by recovery [12], and thus we can assume that δ and δ are in the range [1/4,1/2]. If on the other hand the EVD patient recovers, then it will take longer for patient to be completely clear of the virus. Notice that, based on the ranges given above, the time frames are [6, 16] days from the onset of symptoms of Ebola to death or recovery. The range from the onset of symptoms which commences the course of illness to death was given as 6–16 days in [8], while the range for recovery was cited as 6–11 days [8]. For our baseline parameters, the mean time from the onset of the illness to death or recovery will be in the range of 6–11 days.

Here, we will consider that the recovery rate for quarantined early symptomatic EVD patients lies in the range [0.4829,0.5903] per day with a baseline value of 0.5366 per day, as cited in [16]. Thus r ∈ [0.4829,0.5903]. If we assume that patients quarantined in the hospital have a better chance of surviving than those in the community or at home, without the necessary expert care that some of the quarantined EVD patients may get, then we can consider that the recovery rate for nonquarantined EVD patients, in the community, would be slightly lower. Thus we scale r by some proportion ω ∈ [0,1], so that r = ωr . Since late symptomatic patients have a much lower recovery chance, then both r and r will be lower than r and r , respectively. Hence we will consider that r = κr and r = κr , where 0 < κ ≪ 1. Sometimes late symptomatic EVD patients are removed from the community and quarantined. Here we assume a mean of 2 days so that σ = 0.5 per day.

The fractions θ measure the proportions of individuals moving into various compartments. If we assume that members in the quarantined classes are not left unchecked but have medical professionals checking them and giving them supportive remedies to boost their immune system to fight the Ebola Virus or enable their recovery, then it will be reasonable to assume that θ 6, θ 7, θ 4, and θ 5 are all greater than 0.5. If such an assumption is not made, then the parameters can be chosen to be equal or close to each other.

The parameter Π is chosen to be 555 per day as in [16]. Furthermore, the parameters ρ and ρ and the effective contact rates between probable nonquarantined and, respectively, quarantined individuals and susceptible individuals will be varied to see their effects on the model dynamics. However, the values chosen will be such that the value of R 0 computed is within realistic reported ranges.

The parameter μ, the natural death rate for humans, is chosen based on estimates from [13]. The parameter b measures the time it takes from death to burial of EVD patients. A mean value of 2 days was cited in [14] for the 1995 and 2000 Ebola outbreak epidemics in the Democratic Republic of Congo and Uganda, respectively. For the 2014 West African Outbreak, the estimates were 2.01 days in Liberia and 4.50 days in Sierra Leone [15].

5. Numerical Simulation of the Scaled Reduced Model

The parameter values given in Table 2 were used to carry out some numerical simulations for the reduced model, (81)–(85), when the constant recruitment term is 555 persons per day. The varying parameters were chosen so that R 0 would be within ranges of reported values, which are typically less than 2.5 (see, e.g., [16, 22]). In Figure 2, we show a time series solution for a representative choice of values for the parameters ρ and ρ . Figures 2(a)–2(e) show the long term solutions to the reduced model exhibiting convergence to the stable nontrivial equilibrium in the case where R 0 > 1, the case with sustained infection in the community. Figure 2(f) then shows an example of convergence to the trivial steady state when R 0 < 1, the case where the disease is eradicated. In that example we notice that as R 0 < 1, s → 0 as n → ∞ and this in turn implies that s → 0 as n → ∞ and eventually the system relaxes to the trivial state (s, s , s , r , h) = (1,0, 0,0, 1) for large time.

Figure 2

(a)–(e) Time series showing convergence of the solutions to the steady states for the nondimensional reduced model when the constant recruitment term is 555 persons per day. In this example, θ 1 = 0.85 and ρ = 1.15 and ρ = 0.85, and all the other parameters are as in Table 2, giving values of R 0 = 1.102. In this case, the nonzero steady state is stable and the solution converges to the steady state value as given by (89) as t → ∞. (f) Time series showing the long term behaviour of the variable s in the reduced model for ρ = ρ = 0.8 and all other values of the parameters are as given in Table 2. In this case R 0 = 0.88 and so the only steady state is the trivial steady state which is stable.

We note here that the computed value of R 0 can be shown to be linear in the variables ρ and ρ , when eventually all parameters have been assigned; it may be written as R 0 = r 0 + r 1 ρ + r 2 ρ , where r , i = 0,1, 2, are positive constants that can be shown to be dependent on the other parameters. Thus R 0 will increase linearly with increase in any of the parameters ρ and ρ for fixed given values of the other parameters. Though we have theoretically found, for example, that s becomes infinite when ℛ is near one equivalent to R 0 being near 1/z, this case does not arise because we have assumed in the analysis that z > 1. Thus the case 0 ≤ z < 1 is linked with the trivial steady state.

The situation shown in Figure 3 has important consequences for control strategies. While s varies sharply for a narrow band of reproduction numbers, its values do not change much for larger values of R 0. Referring to Figure 3, an application of a control measure that will reduce R 0 say from a high value of 10 down to 5, a 50% reduction, will not appreciably affect the rest of the disease transmission. Thus the system is best controlled when R 0 is small which can occur in the early stages of the infection or late in the infection when some effective control measures have already been instituted such as effective quarantining or prompt removal of EVD deceased individuals. Notice that as R 0 further increases, the number of susceptible individuals continues to drop. We note, however, that typical values of R 0 computed for the 2014 Ebola outbreak are less than 2.5.

Figure 3

Graph showing the behaviour of the steady states s and s as a function of R 0. (a) This graph shows the form of the steady state s as a function of R 0. (b) This graph shows the form of the steady state s as a function of R 0. The steady state solution s varies greatly only in a narrow band of reproduction numbers but saturates for large values of R 0. On the other hand, s continues to drop to zero as R 0 increases.

Next, we investigate the effect of ξ on R 0 and the model dynamics. In Figure 2(f), we showed an example of convergence to the trivial steady state for ρ = 0.8 = ρ for the nondimensional reduced model when the constant recruitment term is 555 persons per day. The model dynamics yielded an R 0 value of 0.88 < 1, when all other parameters were as given in Table 2. From this scenario, we increased only ξ from its baseline Table 2 value of 0.27 to 0.453. This yields an increase in R 0 to 1.00038 and we see from Figures 4(a) and 4(b) that the disease begins to propagate and stabilize within the community. There is a major peak which starts to decay as EVD deaths begin to rise. An estimated size of the epidemic can be computed as the difference between s and h when the disease dynamics settles to its equilibrium state. As more and more persons become infected, R 0 increases and the estimated size of the epidemic also increases as is expected. In particular, increasing ξ further to 0.48 increases R 0 to 1.01765, and, from Figures 4(c) and 4(d), the difference between s and h is visibly larger compared to the difference from Figures 4(a) and 4(b). Moreover, the oscillatory dynamics becoming more pronounced indicated a higher back and forth movement activity between the s and s and s classes.

Figure 4

(a)–(d) Time series plot showing the propagation and stabilization of EVD to a stable nontrivial steady state for the nondimensional reduced model when the constant recruitment term is 555 persons per day and with ρ = ρ = 0.8 and θ 1 = 0.85 as used in Figure 2(f). Except for ξ that is increased from its baseline value of 0.27, all other parameters are as in Table 2. In graphs (a) and (b), ξ is increased to 0.453. This yields R 0 = 1.00038, slightly bigger than 1. The graphs show that there is a major peak which starts to decay as EVD deaths begin to rise. The size of the epidemic can be estimated as the difference in the areas between the s and h curves as the disease settles to its steady state. Graphs (c) and (d) show the model dynamics when ξ is further increased to 0.48, which yields R 0 = 1.01765. In graphs (c) and (d), the oscillations are more pronounced and the size of the epidemic is larger due to the increased effective contacts with late symptomatic individuals.

In Figures 2(a)–2(e), we showed an example of the long term dynamics of the solutions of the reduced model for ρ = 1.15 and ρ = 0.85. For this case, we obtained R 0 = 1.102 > 1 and the model dynamics show how the reduced model converges to a stable nontrivial equilibrium, when all other parameters are as stated as in Table 2. From this point, if we increase only ξ from its default value of 0.27 to 0.36, we see that R 0 increases to approximately 1.159 and the model exhibits irregular random oscillations with higher frequency but eventually stabilizes (Figures 5(a) and 5(b)). The size of the epidemic is considerably larger in this case. This highlights the importance of reducing contacts between EVD patients and susceptible humans in controlling the size of the disease burden and lowering the impact of the disease.

Figure 5

(a)-(b) Time series plot showing the propagation and stabilization of EVD to a stable nontrivial steady state for the nondimensional reduced model when the constant recruitment term is 555 persons per day and with ρ = 1.15, ρ = 0.85, and θ 1 = 0.85 as in Figures 2(a)–2(e). Except for ξ that is increased from 0.27 to 0.36, all the other parameters are as given in Table 2, and the corresponding R 0 value is R 0 = 1.159. Notice that, in this case, the disease has a higher frequency of oscillations and the difference between the areas under h and s is considerably larger indicating that the size of the disease burden is considerably larger in this case.

5.1. Fade-Outs and Epidemics in Ebola Models

Our model results as highlighted in Figures 2(a)–2(e) , 4, and 5 indicate that it is possible to have a long term endemic situation for Ebola transmission, if conditions are right. In particular, in our model, for the case where we have a relatively large constant recruitment term of 555 persons per day and with available resources to sustain the quarantine efforts then as long as there are people in the community (nonquarantined) with the possibility to come in contact with infectious EVD fluids, then the disease can be sustained as long as R 0 > 1 (Figures 2(a)–2(e), 4, and 5). Increasing control by reducing contacts early enough between suspected and probable individuals with susceptible individuals can bring down the size of R 0 to a value <1 which eventually leads to the eradication of the disease. Thus control which includes quarantining has to be comprehensive and sustained until eradication is achieved.

The 2014 Ebola outbreak did not show sustained disease states. The disease dynamics exhibited epidemic fade-outs. Here, we show that such fade-outs are possible with our model. To investigate the epidemic-like fade-outs, we first note that due to the scaling adopted in our model, our time scales are large. However, any epidemic-like EVD behaviour would be expected to occur over a shorter timescale. Thus, for the results illustrated here, we plot the model dynamics in terms of the original variable by simulating the equations that make up the system, (11)–(23). To illustrate that the large time scales do not affect the long term dynamics of the model results, we first present a graph of the original system in the case where the parameters are maintained as those used in Figure 6 with Π = 555 persons per day, θ 1 = 0.85 and ρ = 1.15 and ρ = 0.85, and with all the other parameters as given in Table 2. The R 0 value was 1.102 and so a sustained disease with no other effort is possible over a long time frame of more than 10,000 days.

Figure 6

(a)–(c) Time series showing convergence of the solutions to the steady states for the full model in dimensional form when the constant recruitment term is 555 persons per day. In this example, θ 1 = 0.85 and ρ = 1.15 and ρ = 0.85, and all the other parameters are as in Table 2, giving a value of R 0 = 1.102. The graph shows the short scale dynamics as well as the long term behaviour showing stability of the nonzero steady state.

When the number of individuals recruited daily reduces, then we can show that, for the case where Π = 3 persons per day, θ 1 = 0.85 and ρ = 2 and ρ = 1, and all the other parameters remain as in Table 2; then an epidemic-like behaviour is obtained (see Figure 7).

Figure 7

(a)-(b) Time series showing epidemic-like behaviours of the solutions to the steady states for the full model in dimensional form when the constant recruitment term is 3 persons per day over a short time scale. In this example, θ 1 = 0.85 and ρ = 2 and ρ = 1, and all the other parameters are as in Table 2. The value of R 0 = 1.63455. The graph shows the short scale dynamics exhibiting epidemic-like behaviour that fades out.

The dynamics of Figures 6 and 7 indicate that quarantining alone is not sufficient to eradicate the Ebola epidemic especially when there is a relative high number of daily recruitment. In fact, quarantining can instead serve as a buffer zone allowing the possibility of sustained disease dynamics when there are a reasonable number of people recruited each day. However, when the daily recruitment is controlled, reduced to a value of 3 per day, then the number of new daily infections is reduced to a low value as depicted in Figures 7(a)–7(c). However, the estimated cumulative number of infections increases daily (see Figure 7(d)).

6. Discussion and Conclusion

In this paper we set out to derive a comprehensive model for the dynamics of Ebola Virus Disease transmission in a complex environment where quarantining is not effective, meaning that some suspected cases escape quarantine while others do not. When the West African countries of Liberia, Guinea, and Sierra Leone came face to face with Ebola Virus Disease infection in 2014, it took the international community some time to react to the crises. As a result, most of the initial cases of EBV infection escaped monitoring and entered the community. African belief systems and other traditional practices further compounded the situation and before long large number of cases of EBV infections were in the community. Even when the international community reacted and started putting in place treatment centres, it still took some time for people to be sensitized on the dangers they are facing. The consequence was that infections continued in families, during funerals, and even in hospitals. People checked into hospitals and would not tell the truth about their case histories and as a result some medical practitioners got exposed to the infection. A case in point is that of Dr. Stella Ameyo Adadevoh, an Ebola victim and everyday hero [31], who prevented the spread of Ebola in Nigeria and paid with her life. We still pay tribute and honour to her and the other health workers whose dedication was inspirational and helpful in curbing the 2014 Ebola outbreak in Africa.

As we now look forward with optimism for a better and Ebola-free tomorrow, there is work going on in the scientific community to develop vaccines [2]. Mathematical modelling of the dynamics and transmission of Ebola provides unique avenues for exploration of possible management scenarios in the event of an EVD outbreak, since, during an outbreak, management of the cases is crucial for containment of the spread of the infection within the community. In this paper we have presented a comprehensive ordinary differential equation model that handles management issues of EVD infection. Our model takes care of quarantine and nonquarantine cases and therefore can be used to predict progression of disease dynamics in the population. Our analysis has shown that the initial response to all suspected cases of EVD infection is crucial. This is captured through the parameter θ 1 which measures the initial fraction of suspected cases that are put into quarantine. We have shown that the basic reproduction number can be indexed by this parameter in the sense that when all cases are initially quarantined the spread of the infection can only take place at the treatment centres, but in cases where all suspected cases escape quarantine, the reproduction number can be large. Our model has been able to quantify the densities of infected and recovered individuals within the population based on baseline parameters identified during the 2014 Ebola Virus Disease outbreak in Africa.

The basic reproduction number in our model depends on the initial exposure rates including exposure to cadavers of EVD victims. The provision of scope for further quarantining during the progression of the infection means that these exposure rates are weighted accordingly, depending on whether or not the system woke up from slumber and picked up those persons who initially escaped quarantine. For example, the parameter τ which measures the effective contact rate between confirmed quarantined early symptomatic individuals and susceptible individuals is eventually scaled by the proportions θ 3 and θ 2 which, respectively, are those proportions of suspected and probable cases that eventually progress to become EVD patients and have escaped quarantine. Thus our framework can progressively be used at each stage to manage the progression of infections in the community.

Our results show that eventually the system settles down to a nonzero fixed point when there is constant recruitment into the population of 555 persons per day and for R 0 > 1. The values of the steady states are completely determined in terms of the parameters in this case. Our analysis also shows that it is possible to control EVD infection in the community provided we reduce and maintain the reproduction number to below unity. Such control measures are possible if there is effective contact tracing and identification of EVD patients and effective quarantining, since a reduction of the proportion of cases that escape quarantine reduces the value of R 0.

Additionally, our model results indicate that when there are a high constant number of recruitment into an EVD community, quarantining alone may not be sufficient to eradicate the disease. It may serve as a buffer enhancing a sustained epidemic. However, reducing the number of persons recruited per day can bring the diseases to very low values.

To demonstrate the feasibility of our results, we performed a pseudoequilibrium approximation to the system derived based on the assumption that the duration of manifestation of EVD infection in the community, per individual, is short when compared with the natural life span of an average human. The reduced model was used to show that all steady state solutions are stable to small perturbations and that there can be oscillatory returns to the equilibrium solution. These results were confirmed with numerical simulations. Given the size of the system, we have not been able to perform a detailed nonlinear analysis on the model. However the discussion on the nature of the parameters for the model is based on the statistics gathered from the 2014 EVD outbreak in Africa and we believe that our model can be useful in characterizing and studying a class of epidemics of Ebola-type. We have not yet carried out a complete sensitivity analysis on all the parameters to determine the most crucial parameters in our model. This is under consideration. Furthermore, the effect of stochasticity seems relevant to study. This and other aspects of the model are under consideration.