Abstract
We propose a new fractional-order model to investigate the transmission and spread of Ebola virus disease. The proposed model incorporates relevant biological factors that characterize Ebola transmission during an outbreak. In particular, we have assumed that susceptible individuals are capable of contracting the infection from a deceased Ebola patient due to traditional beliefs and customs practiced in many African countries where frequent outbreaks of the disease are recorded. We conducted both epidemic and endemic analysis, with a focus on the threshold dynamics characterized by the basic reproduction number. Model parameters were estimated based on the 2014-2015 Ebola outbreak in Sierra Leone. In addition, numerical simulation results are presented to demonstrate the analytical findings.
1. Introduction
In recent decades, fractional calculus theory has been applied in many fields such as mechanical and mechanics, viscoelasticity, bioengineering, finance, optimal theory, optical and thermal system, and electromagnetic field theory [1–3]. Prior studies have shown that fractional calculus is capable of describing rules and development process of some phenomena in natural science [1]. In particular, it has been found that the fractional-order differential system has the advantages of simple modeling, clear parameter meaning, and accurate description for some materials and processes with memory and genetic characteristics [4]. Hence, there is growing interest among researchers to study the role of fractional calculus on modeling real-world problems. One field that has attracted a lot of interest in the application of fractional calculus is mathematical modeling of infectious diseases [2, 3].
In this paper, a fractional-order Ebola epidemic model that incorporates nonlinear incidence rates is proposed and analyzed. A plethora of mathematical models have been proposed to explain, predict as well as quantify the effectiveness of different Ebola virus disease (EVD) intervention strategies since the 2004 when the largest outbreak occurred in Africa (see, for example, [5–14], and references therein). These studies and those cited therein have certainly produced many useful results and improved the existing on Ebola dynamics. One of the limitations of these models, however, is that in most of the studies, the authors were utilizing integer-order mathematical modeling approach except in few recent studies such as [12–14].
In those studies that were based upon fractional calculus, most of the models used the bilinear incidence approach, to describe the spread of the disease. One limitation of the bilinear incidence function is that it assumes that the disease transmission increases whenever the susceptible population increases. This is highly unlikely in practice since an outbreak of any disease is followed by pharmaceutical and nonpharmaceutical interventions. These intervention strategies lead to a saturation in the available population. In mathematical modeling of infectious diseases, the incidence rate is defined as the number of infected individuals per unit time, and it regarded as an important tool for effectively mapping short- and long-term dynamics of the disease [15]. There are several saturated incidence functions in literatures [15, 16]. Among them, the Crowley–Martin function incidence rate has been found to be a very successful tool and has been extensively used [15]. Hence, the framework considered in this paper, will make use of the Crowley–Martin function incidence rate. In addition, for the fractional calculus, we will adopt Caputo’s definition, which is extension of the Riemann–Liouville and has merits on dealing with initial value problem. Another advantage of the Caputo derivative is that its derivative for a constant is zero, while in the Riemann–Liouville sense, the derivative for a constant is nonzero [17].
The paper is organized as follows. In Section 2, a fractional-order Ebola model is formulated and the dynamical analysis is performed. In particular, we conduct both epidemic and endemic analysis, with a focus on the threshold dynamics characterized by the basic reproduction number. Section 3 constitutes numerical results and discussions. Precisely, the proposed model is fitted with the 2014-2015 Ebola data for Sierra Leone. There are also some additional numerical illustrations which have been presented to support analytical findings in the study. Concluding remarks in Section 4 round up the paper.
2. Model Formulation and Analytical Results
2.1. Model Formulation
In this section, we propose a fractional-order Ebola model. The model consists of seven state variables, namely, susceptible individuals , exposed/latently infected individuals , undetected infectious individuals , infectious deceased and undetected Ebola patients , detected and quarantined Ebola patients , deceased Ebola patients in a quarantine facility , and recovered individuals . The recovered population is made up of individuals who have successfully recovered from the infection either naturally or through various health support mechanisms (since the disease has no treatment). These individuals are assumed to recover with long-term immunity which can last for up to 10 years [18]. Hence, the total population in the community at time is given by . The proposed model is governed by the following assumptions:(i)All new recruits into the model are assumed to enter the susceptible class through birth at constant rate ; natural death rate is assumed to occur in all classes at constant rate .(ii)Ebola virus disease (EBOV) is primarily transmitted through direct contact with body fluids, blood of an individual who is sick or has died from Ebola [19]. The force of infection which accounts for disease transmission between susceptible individuals and undetected infectious individuals is given by Here, disease transmission rate is being modeled by the Crowley–Martin incident rate [20]. The parameter is the contact rate due to the infected population, is the inhibition effect due to the susceptible, and is the inhibition effect due to the infected population. One can observe that the function is a saturation function. Thus disease transmission is assumed to be bounded. In particular, the saturation occurs due to disease mitigation strategies carried out by both susceptibles and infectious individuals.(iii)We assume that infectious Ebola patients who succumb to disease or natural-related death while they are not in quarantine facility are capable of generating new infections since their burials are more unlikely to be supervised by health experts. Hence, susceptible individuals are also assumed to acquire infection following contact with a corpse of an infectious deceased and undetected Ebola patient at rate: where is the contact rate between the susceptible and deceased individuals. This assumption is motivated by African practices (the traditional belief systems and customs, for example, washing of deceased individuals) during burial ceremonies.(iv)Upon infection, individuals progress to the latent stage (exposed state) where they incubate the disease for days. Prior studies suggest that this period ranges between 2–21 days. It is worth noting that individuals in this state are not yet infectious; hence, they do not transmit the infection. From the incubation state, infected individuals progress to the undetected infectious state where they are either detected and quarantined at rate day or remain in that state for an average period of days. A proportion of undetected individuals is assumed to successfully recover from infection, and the remainder is assumed to suffer disease-related death. Moreover, a fraction of detected and quarantined individuals is assumed to successfully recover from infection and the complement suffers disease-related death. An individual who succumbs to disease-related death while in quarantine facility is assumed not to be infectious since his/her burial will be supervised by health experts. Further, we assume that it takes an average period of days for a deceased individual to be buried.
Based on the above assumptions, the proposed fractional-order Ebola model is summarized by the following system of equations (Figure 1 illustrates the model structure):with initial conditions:where is the order of the fractional derivative. When , then the model is the integer-order counterpart. The fractional derivative of model (3) is used in the Caputo sense [21], i.e.,where represents the gamma function. The Riemann–Liouville fractional integral of an arbitrary real order of a function is defined by
Model (3) exhibits some problems in what concerns time dimension between left- and right-hand sides of the equations. On the left, the dimension is , whereas on the right-hand side the dimension is . The corrected system corresponding to model (3) is as follows:
From model (7), one can observe that individuals in epidemiological states , , and have no influence on generation of new Ebola cases. Mathematically, we can observe that the variables , , and do not appear in first four equations of (7). Hence, it suffices to investigate the dynamics of the disease based on a reduced system that excludes the last three equations of system (7), i.e.,
2.2. Positivity and Boundedness of Model Solutions
Since model (8) monitors human population, it is imperative to investigate if the model is mathematically and biologically well-poised.
Theorem 1. Let be the unique solution of system (8) for . Then, the solution remains in . Further, the solution is bounded above; that is, where denotes the feasible region and is given bywhere and .
Proof. Based on the results in (10), we conclude that the vector field given by the right-hand side of (8) on each coordinate plane is either tangent to the coordinate plane or points to the interior of . Hence, the domain is a positively invariant region. Moreover, if the initial conditions of system (8) are nonnegative, then it follows that the corresponding solutions of model (8) are nonnegative.
Further, for model (8) to be biological and mathematically meaningful, all model solutions need to be positive and bounded. From (10), we have noted that the model is positively invariant, and in what follows we demonstrate that its solutions are bounded. Since all solutions of model system (8) have been shown to be positively invariant, then it follows that the possible lower-bound for these solutions is zero. Thus, in what follows, we will concentrate on the upper-bound for these solutions.
Let . Adding all the equations of system (8), we haveApplying the Laplace transform leads tothat can be written asApplying the inverse Laplace transform leads towhere . Thus, is bounded from below and above. Hence, one can conclude that the solution is bounded below and above.
2.3. Steady States and Their Stability
In this section, we will present the steady states of model (8) and investigate their stability. In the absence of the disease in the community, that is, , model (8) admits a trivial equilibrium point commonly known as the disease-free equilibrium (DFE), which has and . One of the importance of this trivial equilibrium point is that it is used when computing the model’s basic reproduction number, . The basic reproduction number is an important epidemiological threshold parameter for infectious disease models. It demonstrates the strength of the disease to invade the population. If , it implies that the disease persists. However, if , it implies that disease becomes extinct. By closely following the next-generation matrix method [22], it can be verified that model (8) has a reproduction number as follows:
The threshold quantity, , measures the power of EVD to invade the community. In particular, is the expected number of secondary Ebola cases produced in a completely susceptible population, by one infectious-alive individual during his/her entire infectious period; is defined as the average number of new Ebola infections generated by one infectious corpse introduced in a completely susceptible population.
Theorem 2. (i) For , the disease-free equilibrium of system (8) is globally asymptotically stable whenever . (ii) For , the endemic equilibrium of system (8) is globally asymptotically stable whenever .
Proof. To investigate the global stability of the DFE, we consider the following Lyapunov function:with , , and . Computing the fractional time derivation of along the solutions of system (8), one getsSince at DFE, ; it follows that and if and only if of . Therefore, by Lasalle’s invariance principle [23], the DFE is globally asymptotically stable whenever . This completes the first part of Theorem 2 (i).
To investigate the global stability of the endemic equilibrium, we will make use of Lemma 1.
Lemma 1. (see [24]). Let be continuous and differentiable function with . Then, for any time instant , one has, for all .
Now, we define the following Lyapunov function:where and are constants to be determined. By closely following Lemma 1, we have
At the endemic point, we have the following identities:
Set and .
After some algebraic manipulations, one getswith
Since the arithmetic mean is greater or equal to the geometric mean, it follows that
Therefore, for all and whenever . Further, let , . Note that the function is nonpositive for and if and only if . Now, making use of the properties of the function , we can express and as follows:
Similarly, we have as
Based on the above analysis, we conclude that . Hence, the largest invariant set where is the singleton . Thus, by applying LaSalle’s invariance principle [23], we conclude that the endemic equilibrium point is globally asymptotically stable. This completes the first part of Theorem 2 (ii).
3. Results and Discussion
Our discussion commences with estimation of model parameters and the appropriate fractional order. Thereafter, we will perform some comparison between the fractional-order and integer model. In addition, we will also carry out some numerical illustrations to support some of the analytical findings presented in the study. In order to estimate model parameters for the proposed fractional-order EVD model, the cumulative number of EVD cases reported (sum of confirmed, suspected, and probable cases) in the first 30 weeks of the 2014-2015 Ebola outbreak in Sierra Leone was utilized (see Table 1) [25].
We numerically solved differential equations of system (8) and found the best fractional-order parameter that minimizes the deviations between the real data and model estimates. Precisely, we computed the root-mean-square error (RMSE) to find the optimal parameters:where is the number of observations. Here, is equivalent to 30.
The initial population levels were assumed as follows: , , , and . Model parameters that were drawn from the literature and those that were determined through data fitting are presented in Table 2. Figure 2 illustrates the root-mean-square error (RMSE) for different derivative orders. From the illustration, we can note that the minimum error of estimation for the given data set occurs for . Therefore, according to this illustration, the best fit is obtained by setting . Figure 3(a) shows the model estimates for . In Figure 3(b), the prediction ability between the integer and fractional models is investigated. We can note that the estimates for the integer model are close to the real data for the first 20 weeks there after the estimates significantly deviate from the reported cases compared with the estimates of the fractional-order model. Hence, we conclude that the fractional-order model presents better forecasts compared with the integer model.
(a)
(b)
(a)
(b)
To investigate the effects of different derivative orders, , on the dynamics of the disease, model (8) was simulated with parameter values in Table 2 for , and the results are illustrated in Figure 4. As we can note, with baseline values in Table 2, the disease will persist. In addition, we note that the variables for infected epidemiological classes show that the infection will increase rapidly during first 20 weeks. This is followed by a rapid decline and stability of solutions at the endemic equilibrium. Furthermore, one can observe that as the value of the fractional-order α approaches unity, then the time taken by model variables to converge to their respective equilibrium state increases. In addition, we can also note that for weeks, model solutions for each variable with different fractional-order values converge to a unique solution.
(a)
(b)
(c)
(d)
To support the analytical results in Theorem 2 (i), we simulated model (8) using the following assumed initial conditions and model parameters: , , , , , and (all the other model parameters are as in Table 2, leading to ), and the graphical illustrations are in Figure 5. As one can observe, the illustration supports the analytical results in Theorem 2 (i) that whenever the basic reproduction number is less than unity, then the infection dies out.
(a)
(b)
(c)
(d)
Numerical results in Figures 4 and 5 have shown that the basic reproduction number is an important metric for persistence and extinction of Ebola in the community. Hence, it is imperative that we investigate the influence of model parameters on the growth or decline of . In order to infer on this relationship, we will make use of the approach in [28].
Definition 1. (see [28]). The normalized sensitivity index of , which depends differentiably on a parameter say , is defined by . For example, with , we haveSince the outcome in (28) is positive, it implies that an increase in the magnitude of leads to an increase in the magnitude of . In order to determine the percentage impact, we will use parameter values in Table 2 to compute the numerical values (see Table 3).
Results in Table 3 show that an increase in , , , and will increase the magnitude of . However, an increase in , , , , and will decrease the magnitude of . Most importantly, we note that behavioral change associated with susceptible individuals (modeled by parameter ) has the highest influence on decreasing the magnitude of . Precisely, an increase in by 10% will lead to change in the magnitude of by . In addition, an increase in by 10% will increase by 9.68%.
Results on sensitivity analysis (Table 3) have shown that an increase in case detection (modeled by ) by 10% will reduce the magnitude of by approximately 7.8%. Hence, it imperative that we investigate the threshold detection level that leads to disease extinction. To explore the impact of the case detection on Ebola dynamics, we have varied the values of the case detection parameter , and the results are shown in Figure 6. As we can observe, the higher the detection rate, the less increase in new Ebola cases. In particular, a detection rate of 6.0 per week can sufficiently lead to disease eradication.
(a)
(b)
(c)
4. Concluding Remarks
Mathematical modeling is an important tool for epidemiologists and policy makers since it provides invaluable insights where empirical observations cannot. In this study, we propose a fractional-order modeling framework for Ebola virus disease (EVD) that incorporates nonlinear incidence rates. Recent studies have shown that fractional-order differential equations are more ideal to model many real-world problems in engineering, biology, chemistry, and so on since fractional derivatives are dependent on historical states in addition to the current state and thus they possess memory properties. In addition, the recent mathematical model for infectious disease transmission has shown that the incidence rates are essential functions for tracking and mapping the spread of a disease in a community; hence, any proposed model to track and map the disease during an outbreak needs to incorporate more realistic incidence rates. Majority of mathematical models for EVD in the literature are based on integer-ordinary differential equations which do not have memory properties. Further, a few of those that utilized fractional derivatives disease transmission was modeled using the bilinear incidence which is not realistic since an increase in the number of susceptible population implies that the number of infected individuals per unit time increases. Thus, the primary aim of this study was to develop a more realistic EVD model that incorporates memory effects.
Apart from the inclusion of nonlinear incidences, the model incorporates all relevant biological factors that characterize Ebola transmission and control. In particular, we have assumed that deceased individuals are capable of transmitting the infection. These assumptions were motivated by traditional and cultural belief practiced in most countries where outbreaks of the disease are common. The threshold quantity to determine the number of secondary infections generated by infected individuals (dead or alive) during their entire infectious period was determined. Through sensitivity analysis, we managed to identify model parameters that have a strong influence on either increasing or decreasing the spread of EVD in the community. It was found that parameter which accounts for EVD transmission by infected and alive EVD patients has the greatest influence on increasing the magnitude of the reproduction number. Precisely, an increase in the magnitude of this parameter by 10% may increase the reproduction number by 9.7%. We have also found that disease transmission rate by deceased undetected EVD patients does not have the same impact as . In particular, an increase of by 10% may only increase the reproduction number by 0.3%, and this is far less compared with studies in the literature that utilized bilinear or standard incidence rates. This outcome clearly shows the effects of behavior change by susceptible individuals on the dynamics of EVD during any outbreak.
In addition, by constructing suitable Lyapunov functionals, we have shown that the fractional-order model has a globally asymptotically stable disease-free equilibrium whenever is less than unity. However, if is greater than unity, the fractional-order model has an endemic equilibrium which is globally asymptotically stable. Model parameters were estimated based on the 2014-2015 Ebola outbreak in Sierra Leone. From our numerical results, we noted that for a fractional-order value , the estimates from the proposed model were very close to observed cases than for . Furthermore, we noted that the different values of the fractional-order have no effect on the stability of the endemic equilibrium but influences the time taken for the stability to be attained. In addition, we have also observed that a case detection rate of 6.0 per week could significantly reduce the disease to level close to zero.
This work illustrates the value of mathematical models as a tool to suggest management strategies in disease outbreaks. It is worth noting that the proposed framework is not exhaustive. The model can be extended to incorporate the effects of media campaigns and human mobility.
Data Availability
The data used to support the findings of this study are included within the article.
Conflicts of Interest
The authors declare that they have no conflicts of interest.