# On the Dynamics of a Fractional-Order Ebola Epidemic Model with Nonlinear Incidence Rates

Authors: Mirirai Chinyoka; Tinashe B. Gashirai; Steady Mushayabasa
Journal: Discrete Dynamics in Nature and Society
 (2021)

---

## 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.


---

## Body

## 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 Section2, 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 individualsSt, exposed/latently infected individuals Et, undetected infectious individuals It, infectious deceased and undetected Ebola patients Dt, detected and quarantined Ebola patients Qt, deceased Ebola patients in a quarantine facility Ht, and recovered individuals Rt. The recovered population Rt 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 t is given by Nt=St+Et+It+Dt+Qt+Ht+Rt. 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(1)fS,I=β1SI1+α1S1+α2I.Here, disease transmission rate is being modeled by the Crowley–Martin incident rate [20]. The parameter β1 is the contact rate due to the infected population, α1 is the inhibition effect due to the susceptible, and α2 is the inhibition effect due to the infected population. One can observe that the function fS,I 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:(2)gS,D=β2SD1+α1S,whereβ2 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γ−1 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 −1 or remain in that state for an average period of δ−1 days. A proportion ω of undetected individuals is assumed to successfully recover from infection, and the remainder 1−ω is assumed to suffer disease-related death. Moreover, a fraction f 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 ϕ−1 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 (Figure1 illustrates the model structure):(3)Dtαt0cSt=Λ−fS,I−gS,D−μSDtαt0cEt=fS,I+gS,D−γ+μEDtαt0cIt=γE−μ+θ+δIDtαt0cDt=μ+1−ωδI−ϕDDtαt0cQt=θI−μ+σQDtαt0cHt=1−fσQ−ϕHDtαt0cRt=ωδI+fσQ−μR,with initial conditions:(4)S0=S0E0=E0I0=I0D0=D0Q0=Q0H0=H0R0=R0,where α∈0,1 is the order of the fractional derivative. When α=1, then the model is the integer-order counterpart. The fractional derivative of model (3) is used in the Caputo sense [21], i.e.,(5)Dtαacft=1Γn−α∫0tt−ξn−α−1fnξdξ,n−1<α<n∈ℕ,where Γ represents the gamma function. The Riemann–Liouville fractional integral of an arbitrary real order α>0 of a function ft is defined by(6)Jαft=1Γα∫0tt−ξα−1fξdξ.Figure 1
Model flow diagram.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 time−α, whereas on the right-hand side the dimension is time−1. The corrected system corresponding to model (3) is as follows:(7)Dtαt0cSt=Λα−fS,I−gS,D−μαSDtαt0cEt=fS,I+gS,D−γα+μαEDtαt0cIt=γαE−μα+θα+δαIDtαt0cDt=μα+1−ωδαI−ϕαDDtαt0cQt=θαI−μα+σαQ,Dtαt0cHt=1−fσαQ−ϕαHDtαt0cRt=ωδαI+fσαQ−μαR.From model (7), one can observe that individuals in epidemiological states Qt, Ht, and Rt have no influence on generation of new Ebola cases. Mathematically, we can observe that the variables Qt, Ht, and Rt 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.,(8)Dtαt0cSt=Λα−fS,I−gS,D−μαSDtαt0cEt=fS,I+gS,D−γα+μαEDtαt0cIt=γαE−μα+θα+δαIDtαt0cDt=μα+1−ωδαI−ϕαD.


### 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.
LetXt=St,Et,It,DtT be the unique solution of system (8) for t≥0. Then, the solution Xt remains in ℝ+4. Further, the solution Xt is bounded above; that is, Xt∈Ω where Ω denotes the feasible region and is given by(9)Ω=St,Et,It,Dt∈ℝ+4|0≤Mt≤M˜,where Mt=St+Et+It+Dt+It and M˜=maxM0,Λα/μα.Proof.
(10)Dtαt0cSt|S=0=Λα≥0Dtαt0cEt|E=0=fS,I+gS,D≥0Dtαt0cIt|I=0=γαE≥0Dtαt0cDt|D=0=μα+1−ωδαI≥0.
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 ℝ+4. Hence, the domain ℝ+4 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.
LetMt=St+Et+It+Dt. Adding all the equations of system (8), we have(11)Dtαt0cMt=Dtαt0cSt+Dtαt0cEt+Dtαt0cIt+Dtαt0cDt≤Λα−μαMt.
Applying the Laplace transform leads to(12)sαℒMt−sα−1M0≤Λαs−μαℒMt,that can be written as(13)ℒMt≤Λαs−1sα+μα+M0sα−1sα+μα=Λαsα−1+αsα+μα+M0sα−1sα+μα.
Applying the inverse Laplace transform leads to(14)Mt≤ℒ−1Λαs−1sα+μα+M0sα−1sα+μα+ℒ−1M0sα−1sα+μα≤ΛαtαEα,α+1−μtα+Eα,1−μtαM0≤ΛαμαμαtαEα,α+1−μtα+Eα,1−μtαM0≤maxΛαμα,M0μαtαEα,α+1−μtα+Eα,1−μtα=CΓ1=M˜,where M˜=maxΛα/μα,N˜0. Thus, Mt is bounded from below and above. Hence, one can conclude that the solution Xt 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, Et=It=Dt=0, model (8) admits a trivial equilibrium point commonly known as the disease-free equilibrium (DFE), which has S0=Λα/μα and E0=I0=D0=0. One of the importance of this trivial equilibrium point is that it is used when computing the model’s basic reproduction number, ℛ0. 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 ℛ0>1, it implies that the disease persists. However, if ℛ0<1, 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:(15)ℛ0=ℛ1+ℛ2=β1αγαμα+γαμα+θα+δαΛαμα+α1Λα+β2αγαμα+1−ωαδαϕαμα+γαμα+θα+δαΛαμα+α1Λα.The threshold quantity,ℛ0, measures the power of EVD to invade the community. In particular, ℛ1 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; ℛ2 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α∈0,1, the disease-free equilibrium of system (8) is globally asymptotically stable whenever ℛ0<1. (ii) For α∈0,1, the endemic equilibrium of system (8) is globally asymptotically stable whenever ℛ0>1.Proof.
To investigate the global stability of the DFE, we consider the following Lyapunov function:(16)ℒ0t=β1αγαm1m2+β2αγαm3ϕαm1m2Et+β1αm2+β2αm3ϕαm2It+β2αϕαDt,with m1=μα+γα, m2=μα+θα+δα, and m3=μα+1−ωδα. Computing the fractional time derivation of ℒ0t along the solutions of system (8), one gets(17)Dtαt0cℒ0t=β1αγαm1m2+β2αγαm3ϕαm1m2Dtαt0cEt+β1αm2+β2αm3ϕαm2Dtαt0cIt+β2αϕαDtαt0cDt=β1αγαm1m2+β2αγαm3ϕαm1m2β1αIS1+α1S1+α2I+β2αDS1+α1S−β1αI+β2αD≤β1αγαm1m2+β2αγαm3ϕαm1m2β1αI+β2αDS1+α1S−β1αI+β2αD≤β1αγαm1m2Λαμα+α1Λα+β2αγαm3ϕαm1m2Λαμα+α1Λα−1β1αI+β2αD=ℛ0−1β1αI+β2αD.
Since at DFE,ℛ0<1; it follows that Dtαt0cℒ0t<0 and Dtαt0cℒ0t=0 if and only if ℛ0=1 of It=Dt=0. Therefore, by Lasalle’s invariance principle [23], the DFE is globally asymptotically stable whenever ℛ0<1. This completes the first part of Theorem 2 (i).
To investigate the global stability of the endemic equilibrium, we will make use of Lemma1.Lemma 1.
(see [24]). Let xt be continuous and differentiable function with xt∈ℛ+. Then, for any time instant t≥t0, one has(18)Dtαxt−x∗−x∗lnxtx∗≤1−x∗xtDtαxt,x∗∈ℛ+, for all α∈0,1.Now, we define the following Lyapunov function:(19)ℒ1t=S−S∗−S∗lnSS∗+E−E∗−E∗lnEE∗+a1I−I∗−I∗lnII∗+a2D−D∗−D∗lnDD∗,where a1 and a2 are constants to be determined. By closely following Lemma 1, we have(20)Dtαt0cℒ1t≤1−S∗StDtαt0cSt+1−E∗EtDtαt0cEt+a11−I∗ItDtαt0cIt+a21−D∗DtDtαt0cDt.At the endemic point, we have the following identities:(21)Λα=fS∗,I∗+gS∗,D∗+μαS∗,m1=fS∗,I∗E∗+gS∗,D∗E∗,m2=γαE∗I∗,m3=ϕαD∗I∗.Seta1=fS∗,I∗+gS∗,D∗/γαE∗ and a2=fS∗,I∗+gS∗,D∗/ϕαD∗.After some algebraic manipulations, one gets(22)Dtαℒt≤μS∗W1+fS∗,I∗W2+gS∗,D∗W3,with(23)W1=2−S∗St−StS∗,W2=4+S∗fSt,ItStfS∗,I∗−S∗St−E∗fSt,ItEtfS∗,I∗−EtI∗E∗It−DtD∗−D∗ItDtI∗,W3=4+S∗gSt,DtStgS∗,D∗−S∗St−E∗gSt,DtEtgS∗,D∗−EtI∗E∗It−DtD∗−D∗ItDtI∗.Since the arithmetic mean is greater or equal to the geometric mean, it follows that(24)S∗St+StS∗≥2S∗St⋅StS∗.Therefore,W1<0 for all St≠S∗ and W1=0 whenever St=S∗. Further, let Φx=1−x+lnx, x>0. Note that the function Φx is nonpositive for x>0 and Φx=0 if and only if x=0. Now, making use of the properties of the function Φx, we can express W2 and W3 as follows:(25)W2=4+S∗fSt,ItStfS∗,I∗−S∗St−E∗fSt,ItEtfS∗,I∗−EtI∗E∗It−DtD∗−D∗ItDtI∗=ΦS∗St+ΦE∗fSt,ItEtfS∗,I∗+ΦEtI∗E∗It+ΦD∗ItDtI∗+S∗fSt,ItStfS∗,I∗−DtD∗−lnS∗fSt,ItD∗StfS∗,I∗Dt≤S∗fSt,ItStfS∗,I∗−DtD∗−lnS∗fSt,ItD∗StfS∗,I∗Dt=S∗fSt,ItStfS∗,I∗−lnS∗fSt,ItStfS∗,I∗−DtD∗+lnDtD∗≤0.Similarly, we haveW3 as(26)W3=4+S∗gSt,DtStgS∗,D∗−S∗St−E∗gSt,DtEtgS∗,D∗−EtI∗E∗It−DtD∗−D∗ItDtI∗=ΦS∗St+ΦE∗fSt,DtEtfS∗,D∗+ΦEtI∗E∗It+ΦD∗ItDtI∗+S∗fSt,DtStfS∗,D∗−DtD∗−lnS∗D∗fSt,DtStDtfS∗,D∗≤−lnS∗D∗fSt,DtStDtfS∗,D∗+S∗fSt,DtStfS∗,D∗−DtD∗≤S∗fSt,DtStfS∗,D∗−lnS∗fSt,DtStfS∗,D∗−lnD∗Dt−DtD∗=S∗fSt,DtStfS∗,D∗−lnS∗fSt,DtStfS∗,D∗−DtD∗+lnDtD∗≤0.Based on the above analysis, we conclude thatDtαℒ1t≤0. Hence, the largest invariant set where Dtαℒ1t=0 is the singleton E∗. 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).


## 2.1. Model Formulation

In this section, we propose a fractional-order Ebola model. The model consists of seven state variables, namely, susceptible individualsSt, exposed/latently infected individuals Et, undetected infectious individuals It, infectious deceased and undetected Ebola patients Dt, detected and quarantined Ebola patients Qt, deceased Ebola patients in a quarantine facility Ht, and recovered individuals Rt. The recovered population Rt 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 t is given by Nt=St+Et+It+Dt+Qt+Ht+Rt. 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(1)fS,I=β1SI1+α1S1+α2I.Here, disease transmission rate is being modeled by the Crowley–Martin incident rate [20]. The parameter β1 is the contact rate due to the infected population, α1 is the inhibition effect due to the susceptible, and α2 is the inhibition effect due to the infected population. One can observe that the function fS,I 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:(2)gS,D=β2SD1+α1S,whereβ2 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γ−1 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 −1 or remain in that state for an average period of δ−1 days. A proportion ω of undetected individuals is assumed to successfully recover from infection, and the remainder 1−ω is assumed to suffer disease-related death. Moreover, a fraction f 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 ϕ−1 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 (Figure1 illustrates the model structure):(3)Dtαt0cSt=Λ−fS,I−gS,D−μSDtαt0cEt=fS,I+gS,D−γ+μEDtαt0cIt=γE−μ+θ+δIDtαt0cDt=μ+1−ωδI−ϕDDtαt0cQt=θI−μ+σQDtαt0cHt=1−fσQ−ϕHDtαt0cRt=ωδI+fσQ−μR,with initial conditions:(4)S0=S0E0=E0I0=I0D0=D0Q0=Q0H0=H0R0=R0,where α∈0,1 is the order of the fractional derivative. When α=1, then the model is the integer-order counterpart. The fractional derivative of model (3) is used in the Caputo sense [21], i.e.,(5)Dtαacft=1Γn−α∫0tt−ξn−α−1fnξdξ,n−1<α<n∈ℕ,where Γ represents the gamma function. The Riemann–Liouville fractional integral of an arbitrary real order α>0 of a function ft is defined by(6)Jαft=1Γα∫0tt−ξα−1fξdξ.Figure 1
Model flow diagram.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 time−α, whereas on the right-hand side the dimension is time−1. The corrected system corresponding to model (3) is as follows:(7)Dtαt0cSt=Λα−fS,I−gS,D−μαSDtαt0cEt=fS,I+gS,D−γα+μαEDtαt0cIt=γαE−μα+θα+δαIDtαt0cDt=μα+1−ωδαI−ϕαDDtαt0cQt=θαI−μα+σαQ,Dtαt0cHt=1−fσαQ−ϕαHDtαt0cRt=ωδαI+fσαQ−μαR.From model (7), one can observe that individuals in epidemiological states Qt, Ht, and Rt have no influence on generation of new Ebola cases. Mathematically, we can observe that the variables Qt, Ht, and Rt 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.,(8)Dtαt0cSt=Λα−fS,I−gS,D−μαSDtαt0cEt=fS,I+gS,D−γα+μαEDtαt0cIt=γαE−μα+θα+δαIDtαt0cDt=μα+1−ωδαI−ϕαD.


## 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.
LetXt=St,Et,It,DtT be the unique solution of system (8) for t≥0. Then, the solution Xt remains in ℝ+4. Further, the solution Xt is bounded above; that is, Xt∈Ω where Ω denotes the feasible region and is given by(9)Ω=St,Et,It,Dt∈ℝ+4|0≤Mt≤M˜,where Mt=St+Et+It+Dt+It and M˜=maxM0,Λα/μα.Proof.
(10)Dtαt0cSt|S=0=Λα≥0Dtαt0cEt|E=0=fS,I+gS,D≥0Dtαt0cIt|I=0=γαE≥0Dtαt0cDt|D=0=μα+1−ωδαI≥0.
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 ℝ+4. Hence, the domain ℝ+4 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.
LetMt=St+Et+It+Dt. Adding all the equations of system (8), we have(11)Dtαt0cMt=Dtαt0cSt+Dtαt0cEt+Dtαt0cIt+Dtαt0cDt≤Λα−μαMt.
Applying the Laplace transform leads to(12)sαℒMt−sα−1M0≤Λαs−μαℒMt,that can be written as(13)ℒMt≤Λαs−1sα+μα+M0sα−1sα+μα=Λαsα−1+αsα+μα+M0sα−1sα+μα.
Applying the inverse Laplace transform leads to(14)Mt≤ℒ−1Λαs−1sα+μα+M0sα−1sα+μα+ℒ−1M0sα−1sα+μα≤ΛαtαEα,α+1−μtα+Eα,1−μtαM0≤ΛαμαμαtαEα,α+1−μtα+Eα,1−μtαM0≤maxΛαμα,M0μαtαEα,α+1−μtα+Eα,1−μtα=CΓ1=M˜,where M˜=maxΛα/μα,N˜0. Thus, Mt is bounded from below and above. Hence, one can conclude that the solution Xt 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, Et=It=Dt=0, model (8) admits a trivial equilibrium point commonly known as the disease-free equilibrium (DFE), which has S0=Λα/μα and E0=I0=D0=0. One of the importance of this trivial equilibrium point is that it is used when computing the model’s basic reproduction number, ℛ0. 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 ℛ0>1, it implies that the disease persists. However, if ℛ0<1, 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:(15)ℛ0=ℛ1+ℛ2=β1αγαμα+γαμα+θα+δαΛαμα+α1Λα+β2αγαμα+1−ωαδαϕαμα+γαμα+θα+δαΛαμα+α1Λα.The threshold quantity,ℛ0, measures the power of EVD to invade the community. In particular, ℛ1 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; ℛ2 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α∈0,1, the disease-free equilibrium of system (8) is globally asymptotically stable whenever ℛ0<1. (ii) For α∈0,1, the endemic equilibrium of system (8) is globally asymptotically stable whenever ℛ0>1.Proof.
To investigate the global stability of the DFE, we consider the following Lyapunov function:(16)ℒ0t=β1αγαm1m2+β2αγαm3ϕαm1m2Et+β1αm2+β2αm3ϕαm2It+β2αϕαDt,with m1=μα+γα, m2=μα+θα+δα, and m3=μα+1−ωδα. Computing the fractional time derivation of ℒ0t along the solutions of system (8), one gets(17)Dtαt0cℒ0t=β1αγαm1m2+β2αγαm3ϕαm1m2Dtαt0cEt+β1αm2+β2αm3ϕαm2Dtαt0cIt+β2αϕαDtαt0cDt=β1αγαm1m2+β2αγαm3ϕαm1m2β1αIS1+α1S1+α2I+β2αDS1+α1S−β1αI+β2αD≤β1αγαm1m2+β2αγαm3ϕαm1m2β1αI+β2αDS1+α1S−β1αI+β2αD≤β1αγαm1m2Λαμα+α1Λα+β2αγαm3ϕαm1m2Λαμα+α1Λα−1β1αI+β2αD=ℛ0−1β1αI+β2αD.
Since at DFE,ℛ0<1; it follows that Dtαt0cℒ0t<0 and Dtαt0cℒ0t=0 if and only if ℛ0=1 of It=Dt=0. Therefore, by Lasalle’s invariance principle [23], the DFE is globally asymptotically stable whenever ℛ0<1. This completes the first part of Theorem 2 (i).
To investigate the global stability of the endemic equilibrium, we will make use of Lemma1.Lemma 1.
(see [24]). Let xt be continuous and differentiable function with xt∈ℛ+. Then, for any time instant t≥t0, one has(18)Dtαxt−x∗−x∗lnxtx∗≤1−x∗xtDtαxt,x∗∈ℛ+, for all α∈0,1.Now, we define the following Lyapunov function:(19)ℒ1t=S−S∗−S∗lnSS∗+E−E∗−E∗lnEE∗+a1I−I∗−I∗lnII∗+a2D−D∗−D∗lnDD∗,where a1 and a2 are constants to be determined. By closely following Lemma 1, we have(20)Dtαt0cℒ1t≤1−S∗StDtαt0cSt+1−E∗EtDtαt0cEt+a11−I∗ItDtαt0cIt+a21−D∗DtDtαt0cDt.At the endemic point, we have the following identities:(21)Λα=fS∗,I∗+gS∗,D∗+μαS∗,m1=fS∗,I∗E∗+gS∗,D∗E∗,m2=γαE∗I∗,m3=ϕαD∗I∗.Seta1=fS∗,I∗+gS∗,D∗/γαE∗ and a2=fS∗,I∗+gS∗,D∗/ϕαD∗.After some algebraic manipulations, one gets(22)Dtαℒt≤μS∗W1+fS∗,I∗W2+gS∗,D∗W3,with(23)W1=2−S∗St−StS∗,W2=4+S∗fSt,ItStfS∗,I∗−S∗St−E∗fSt,ItEtfS∗,I∗−EtI∗E∗It−DtD∗−D∗ItDtI∗,W3=4+S∗gSt,DtStgS∗,D∗−S∗St−E∗gSt,DtEtgS∗,D∗−EtI∗E∗It−DtD∗−D∗ItDtI∗.Since the arithmetic mean is greater or equal to the geometric mean, it follows that(24)S∗St+StS∗≥2S∗St⋅StS∗.Therefore,W1<0 for all St≠S∗ and W1=0 whenever St=S∗. Further, let Φx=1−x+lnx, x>0. Note that the function Φx is nonpositive for x>0 and Φx=0 if and only if x=0. Now, making use of the properties of the function Φx, we can express W2 and W3 as follows:(25)W2=4+S∗fSt,ItStfS∗,I∗−S∗St−E∗fSt,ItEtfS∗,I∗−EtI∗E∗It−DtD∗−D∗ItDtI∗=ΦS∗St+ΦE∗fSt,ItEtfS∗,I∗+ΦEtI∗E∗It+ΦD∗ItDtI∗+S∗fSt,ItStfS∗,I∗−DtD∗−lnS∗fSt,ItD∗StfS∗,I∗Dt≤S∗fSt,ItStfS∗,I∗−DtD∗−lnS∗fSt,ItD∗StfS∗,I∗Dt=S∗fSt,ItStfS∗,I∗−lnS∗fSt,ItStfS∗,I∗−DtD∗+lnDtD∗≤0.Similarly, we haveW3 as(26)W3=4+S∗gSt,DtStgS∗,D∗−S∗St−E∗gSt,DtEtgS∗,D∗−EtI∗E∗It−DtD∗−D∗ItDtI∗=ΦS∗St+ΦE∗fSt,DtEtfS∗,D∗+ΦEtI∗E∗It+ΦD∗ItDtI∗+S∗fSt,DtStfS∗,D∗−DtD∗−lnS∗D∗fSt,DtStDtfS∗,D∗≤−lnS∗D∗fSt,DtStDtfS∗,D∗+S∗fSt,DtStfS∗,D∗−DtD∗≤S∗fSt,DtStfS∗,D∗−lnS∗fSt,DtStfS∗,D∗−lnD∗Dt−DtD∗=S∗fSt,DtStfS∗,D∗−lnS∗fSt,DtStfS∗,D∗−DtD∗+lnDtD∗≤0.Based on the above analysis, we conclude thatDtαℒ1t≤0. Hence, the largest invariant set where Dtαℒ1t=0 is the singleton E∗. 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 Table1) [25].Table 1
Cumulative number of Ebola cases for Sierra Leone for the duration 2014-2015.
Week123456789101112131415Cases1111223344566612Week161718192021222324252627282930Cases1212121212508189136158252337442525574We 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:(27)RMSE=∑i=1nreal data−estimated cases2n,where n is the number of observations. Here, n is equivalent to 30.The initial population levels were assumed as follows:S0=6900, E0=0, I0=1, and D0=0. 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 α=0.937635. Therefore, according to this illustration, the best fit is obtained by setting α=0.937635. Figure 3(a) shows the model estimates for α=0.937635. 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.Table 2
Description of model parameters of system (8), indicating baseline, ranges, and references.
SymbolDescriptionΛRecruitment rate11826 week−1[26]μNatural death rate0.00054 week−1[26]1/γIncubation period1.498 weeks[26]θHospitalization rate1.58 week−1[27]1/δTime from infection to death2.23 weeks[27]1/ϕMean time from death to Ebola to burial0.762 weeks[26]1−ωCase fatality ratio for unhospitalized patients0.75[27]β1Infection rate through direct contact with a nondeceased infectious patient8.25×10−5Fittingβ2Infection rate through direct contact with a deceased Ebola patient10−5Fittingα1Inhibition effect due to the susceptible individuals9.5×10−6Fittingα2Inhibition effect due to the infectious individuals2.5×10−9FittingFigure 2
The root-mean-square error (RMSE) for different derivative orders. In (a) a more wider range ofα is considered and in (b) we refined the range of α so as to improve precision on determining the minimum error of estimation and we observed that it occurs when α=0.937635.
(a)(b)Figure 3
(a) The estimation of the fractional-order model withα=0.937635. (b) Prediction of the fractional-order (with α=0.937635) and integer model.
(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 α=0.8,0.9,093,1.0, 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 t≥30 weeks, model solutions for each variable with different fractional-order values converge to a unique solution.Figure 4
Dynamical solutions of model (8) for different fractional-order values. The initial population levels were assumed as follows: S0=6900, E0=0, I0=1, and D0=0. All model parameters are as in Table 2, and this yields ℛ0=4.3857.
(a)(b)(c)(d)To support the analytical results in Theorem2 (i), we simulated model (8) using the following assumed initial conditions and model parameters: S0=6900, E0=100, I0=50, D0=0, β1=8.25×10−6, and β2=4×10−5 (all the other model parameters are as in Table 2, leading to ℛ0=0.9564), 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.Figure 5
Dynamical solutions of model (8) for different fractional-order values. The initial population levels were assumed as follows: S0=6900, E0=100, I0=50, and D0=0. We set β1=8.25×10−6 and β2=4×10−5, and all the other model parameters are as in Table 2, leading to ℛ0=0.9564.
(a)(b)(c)(d)Numerical results in Figures4 and 5 have shown that the basic reproduction number ℛ0 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 ℛ0. 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 ℛ0, which depends differentiably on a parameter say r, is defined by Γrℛ0=∂ℛ0/∂r×r/ℛ0. For example, with r=Λ, we have(28)ΓΛℛ0≔∂ℛ0∂Λ×Λℛ0=μμ+α1Λ.
Since the outcome in (28) is positive, it implies that an increase in the magnitude of Λ leads to an increase in the magnitude of ℛ0. In order to determine the percentage impact, we will use parameter values in Table 2 to compute the numerical values (see Table 3).Table 3
Sensitivity index.
Parameterα1β1θβ2ϕωμΛγSensitivity index−0.9950.97−0.780.03−0.03−0.01−0.0060.0050.0008Results in Table3 show that an increase in β1, β2, Λ, and γ will increase the magnitude of ℛ0. However, an increase in α1, θ, ϕ, ω, and μ will decrease the magnitude of ℛ0. Most importantly, we note that behavioral change associated with susceptible individuals (modeled by parameter α1) has the highest influence on decreasing the magnitude of ℛ0. Precisely, an increase in α1 by 10% will lead to change in the magnitude of ℛ0 by 9.95%. In addition, an increase in β1 by 10% will increase ℛ0 by 9.68%.Results on sensitivity analysis (Table3) have shown that an increase in case detection (modeled by θ) by 10% will reduce the magnitude of ℛ0 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.Figure 6
Effects of case detection on Ebola dynamics.
(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β1 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 β2 does not have the same impact as β1. In particular, an increase of β2 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ℛ0 is less than unity. However, if ℛ0 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 α=0.937635, the estimates from the proposed model were very close to observed cases than for α=1. 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.



---

*Source: 2125061-2021-12-03.xml*