Computer Modeling in Engineering & Sciences

DOI: 10.32604/cmes.2021.015472

ARTICLE

Modeling Dysentery Diarrhea Using Statistical Period Prevalence

Fouad A. Abolaban*

Department of Nuclear Engineering, College of Engineering, King Abdulaziz University, Jeddah, 21589, Saudi Arabia
*Corresponding Author: Fouad A. Abolaban. Email: fabolaban@kau.edu.sa
Received: 21 December 2020; Accepted: 07 February 2021

1  Introduction

The literature provides mathematical models for the transmission of infectious diseases. These models play a significant role in quantifying and evaluating the effective control and preventive measures of infectious diseases [1–3]. Furthermore, mathematical modeling has been used in several ways as a versatile and effective way of studying the dynamics of infectious disease transmission. This can include the classic susceptible, infected, and recovered (SIR) model or more advanced models [4]. Mathematical analysis and numerical simulations can be collectively used for the development and evaluation of persuasive control measures.

It is popularly known that mathematical models can predict the emergence of infectious diseases and epidemics, which are beneficial for public health planning and initiatives. By using compartmental models as a simple mathematical structure, the complex dynamics of epidemiological processes can be examined [5]. These compartmental models divide the population into two distinct health categories. The first category is depicted by S, and involves those susceptible to pathogen infection, while the second, denoted by I, involves the pathogen-infected humans. The manner in which these two populations interact is based on phenomenological assumptions used to develop the models. Ordinary differential equations (ODEs) are typically used to develop these models.

Additionally, other populations denoted by R, which is the image of the immune/removed/ recovered compartment, are considered to make these models more practical. A significant challenge here is to obtain sufficient parameters for a particular disease, which would determine the factors affecting potential control measures, such as medication or vaccination. The crucial question is about the execution of such measures from an optimal viewpoint. Several notable attempts have recently been made to introduce this research program for various diseases with integer compartmental models [6–11].

Over the past few decades, many scientists have shown that fractional models can effectively represent natural phenomena compared to integer-order differential equations. Therefore, fractional calculus has gained more importance and popularity for modeling realistic cases, especially memory effects [12–14]. Due to this particular function, various fractional operators have been developed to accurately model the memory effects of various types of diseases [15–22]. Nevertheless, further research is required to explain such complex dynamics. Classical fractional models with singular operators cannot effectively model the non-locality of real-world applications. To overcome this challenge, we investigated and examined a new fractional version of an epidemiological model for the dysentery diarrhea involving the ABC operator, known to have a non-singular kernel with memory effects.

2  Formulation of the Dysentery Diarrhea System

This section demonstrates the formulation of the dysentery diarrhea disease in terms of a deterministic model based on a nonlinear system of ODEs over a finite time interval [0, T], 0 < t < T. The model comprises four population classes. Three classes are reserved for samples from the human population, including those that can be infected with the disease called susceptible S(t). In addition, those that transmit the disease are denoted by I(t), and those that successfully recover from the disease are denoted by R(t). The other class is the dynamics for the pathogen population (concentration of Shigella dysenteriae), denoted by B(t).

The model is designed based on the following assumptions:

•   Transmission of the dysentery diarrhea disease occurs through multiple pathways.

•   There is a homogeneously mixed population

•   Π is the rate of recruitment of those susceptible either by birth or immigration

•   Standard incidence is assumed in the human to human interaction

•   Logistic phenomenon is taken into consideration in the human to environment interaction, which is represented below:

Πh=βhINΠB=βeBB+C,(1)

where C stands for the Shigella concentration that causes a 25% likelihood of getting the disease

•   βh and βe stand for rates of ingesting Shigella from human to human interaction and through a contaminated environment, respectively. βh=ap (a = contact rate and p = probability of disease transmission per contact)

•   After losing immunity, individuals return to S(t) at a rate of α

•   Infected individuals cause concentration of Shigella at a rate of ε

•   Shigella population dies at a rate of σ

•   Rate of recovery of infected humans is γ

•   Natural mortality rate of humans is μ where death due to dysentery diarrhea disease occurs at a rate of d

•   Non-negativity is assumed for all biological parameters introduced within the model

Hence, after incorporating all the above assumptions and considering the Atangana–Baleanu differential operator taken in the Caputo sense [15], we have obtained the following coupled nonlinear system of ordinary differential equations:

ABC(ν)1-ν∫0tS′(ζ)Eν[-ν1-ν(t-ζ)ν]dζ=Πν+ανR(t)-(Πh+ΠB+μν)S(t),ABC(ν)1-ν∫0tI′(ζ)Eν[-ν1-ν(t-ζ)ν]dζ=(Πh+ΠB)S(t)-(μν+γν+dν)I(t),ABC(ν)1-ν∫0tR′(ζ)Eν[-ν1-ν(t-ζ)ν]dζ=γνI(t)-(μν+αν)R(t),ABC(ν)1-ν∫0tB′(ζ)Eν[-ν1-ν(t-ζ)ν]dζ=ενI(t)-σνB(t),(2)

subject to the following initial conditions:

where Πh=βhνI(t)N(t),ΠB=βeνB(t)B(t)+C, and N(t) = S(t)+I(t)+R(t).

3  Existence and Uniqueness

This section presents the existence and uniqueness of solutions of the proposed model using the techniques of fixed point theory. Here, we denote E=C([0,T],ℝ), the Banach space of all continuous real-valued function equipped with the norm defined by:

where

Thus, the proposed fractional model takes the forms shown below:

S(t)-S(0)=ABCI0+ν{Πν+ανR(t)-(Πh+ΠB+μν)S(t))},I(t)-I(0)=ABCI0+ν{(Πh+ΠB)S(t)-(μν+γν+dν)I(t)},R(t)-R(0)=ABCI0+ν{γνI(t)-(μν+αν)R(t)},B(t)-B(0)=ABCI0+ν{ενI(t)-σνB(t)}.(3)

We then obtain:

S(t)=S(0)+1-νN(ν)F1(t,S(t))+νN(ν)1Γ(ν)∫0t(t-x)ν-1F1(x,S(x))dx,I(t)=I(0)+1-νN(ν)F2(t,I(t))+νN(ν)1Γ(ν)∫0t(t-x)ν-1F2(x,I(x))dx,R(t)=R(0)+1-νN(ν)F3(t,R(t))+νN(ν)1Γ(ν)∫0t(t-x)ν-1F3(x,R(x))dx,B(t)=B(0)+1-νN(ν)F4(t,B(t))+νN(ν)1Γ(ν)∫0t(t-x)ν-1F4(x,B(x))dx,(4)

where,

F1(t,S(t))=Πν+ανR(t)-(Πh+ΠB+μν)S(t),F2(t,I(t))=(Πh+ΠB)S(t)-(μν+γν+dν)I(t),F3(t,R(t))=γνI(t)-(μν+αν)R(t),F4(t,B(t))=ενI(t)-σνB(t).(5)

The kernels in Eq. (5) satisfy the Lipschitz condition for 0≤Mi<0, i=1,2,⋯4. When S(t) and S*(t) are two functions, we get:

∥F1(t,S(t))-F1(t,S*(t))∥=∥Πν+ανR(t)-(Πh+ΠB+μν)S(t)-(Πν+ανR(t)-(Πh+ΠB+μν)S*(t))∥=∥(Πh+ΠB+μν)(S*(t)-S(t))∥≤(Πh+ΠB+μν)∥S(t)-S*(t)∥=M1∥S(t)-S*(t)∥,(6)

where M1=(Πh+ΠB+μν).

Thus,

∥F1(t,S(t))-F1(t,S*(t))∥≤M1∥S(t)-S*(t)∥,(7)

Repeating the same procedure above, yields:

∥F2(t,E(t))-F2(t,E*(t))∥≤M2∥E(t)-E*(t)∥,∥F3(t,Q(t))-F3(t,Q*(t))∥≤M3∥Q(t)-Q*(t)∥,∥F4(t,IA(t))-F4(t,IA*(t))∥≤M4∥IA(t)-IA*(t)∥.(8)

Subsequently, Eq. (4) gives:

Sn(t)=S(0)+1-νN(ν)F1(t,Sn-1(t))+νN(ν)1Γ(ν)∫0t(t-x)ν-1F1(x,Sn-1(x))dx,In(t)=E(0)+1-νN(ν)F2(t,In-1(t))+νN(ν)1Γ(ν)∫0t(t-x)ν-1F2(x,In-1(x))dx,Rn(t)=Q(0)+1-νN(ν)F3(t,Rn-1(t))+νN(ν)1Γ(ν)∫0t(t-x)ν-1F3(x,Rn-1(x))dx,Bn(t)=B(0)+1-νN(ν)F4(t,Bn-1(t))+νN(ν)1Γ(ν)∫0t(t-x)ν-1F4(x,Bn-1(x))dx,(9)

where S(t)≥0, I(t)≥0, R(t)≥0, B(t)≥0. The difference between successive components can be denoted by Φni,i=1,2,⋯6, respectively. Thus in view of Eqs. (7)–(9), we obtain:

∥Φn1(t)∥=1-νN(ν)M1∥S(t)-S*(t)∥+νN(ν)M1Γ(ν)∫0t(t-x)ν-1∥S(t)-S*(t)∥dx,∥Φn2(t)∥=1-νN(ν)M2∥I(t)-I*(t)∥+νN(ν)M2Γ(ν)∫0t(t-x)ν-1∥I(t)-I*(t)∥dx,∥Φn3(t)∥=1-νN(ν)M3∥R(t)-R*(t)∥+νN(ν)M3Γ(ν)∫0t(t-x)ν-1∥R(t)-R*(t)∥dx,∥Φn4(t)∥=1-νN(ν)M4∥B(t)-B*(t)∥+νN(ν)M4Γ(ν)∫0t(t-x)ν-1∥B(t)-B*(t)∥dx,(10)

Theorem 3.1 The fractional proposed model possesses a unique solution for t∈[0,T] if the condition is satisfied

(1-νN(ν)Mi+1N(ν)Miγ(ν)Tν)<1,i=1,2,…,6.(11)

Proof. Based on the assumptions that S(t), E(t), Q(t), IA(t), IS(t), R(t) are bounded functions, it is therefore clear that the kernels F1, F2, F3, F4, F5, F6 from Eqs. (7)–(8) satisfy the Lipschitz condition. Hence, Eq. (10) can be viewed as:

∥Φn1(t)∥≤(1-νN(ν)M1+1N(ν)M1Γ(ν)Tν)n,∥Φn2(t)∥≤(1-νN(ν)M2+1N(ν)M2Γ(ν)Tν)n,∥Φn3(t)∥≤(1-νN(ν)M3+1N(ν)M3Γ(ν)Tν)n,∥Φn4(t)∥≤(1-νN(ν)M4+1N(ν)M4Γ(ν)Tν)n.(12)

Hence, the sequences above exist as n→∞,∥Φni(t)∥→0,i=1,2,⋯6. Also, using the triangular inequality for any k value, Eq. (12) yields:

∥Sn+k(t)-Sn(t)∥≤∑i=n+1n+kP1i=P1n+1-P1n+k+11-P1,∥In+k(t)-In(t)∥≤∑i=n+1n+kP2i=P2n+1-P2n+k+11-P2,∥Rn+k(t)-Rn(t)∥≤∑i=n+1n+kP3i=P3n+1-P3n+k+11-P3,∥Bn+k(t)-Bn(t)∥≤∑i=n+1n+kP4i=P4n+1-P4n+k+11-P4,(13)

4  Stability and Basic Reproductive Number

In this section, the stability of the ailment-free equilibrium and its analytic conditions will be discussed. From the derivation of the existence of equilibria in [23], the basic reproduction number in the fractional model is given by:

R0=βhνμν+γν+dν+(Πβeε)νμνσν(μν+γν+dν)C.(14)

Theorem 4.1 The ailment-free equilibrium E0 is locally asymptotically stable if R0≤1 and unstable if R0>1.

Proof. Through the concept of the Jacobian matrix, local stability at E0 can be achieved by:

J[E0]=[-μν-βhνφν-βBΠKμ0βhν-(μν+dν+γν)0βBΠKμ0γν-(φν+μ)00ε0-σν.](15)

The associated eigenvalues are Π1=-μν,Π2=-(μν+φν) and the solutions of Π2+((μν+dν+γν+σν)-βh)Π+(μν+dν+γν)σν(1-R0)=0. Now, if R0<1, ((μν+dν+γν+σν)-βh)<0 and ((μν+dν+γν+σν)-βh)(μν+dν+γν)σν(1-R0)>0, in accordance with the Hurwitz principle, the expression Π2+((μν+dν+γν+σν)-βh)Π+(μν+dν+γν)σν(1-R0)=0, possesses negative real eigenvalues. Thus, E0 is locally asymptotically stable. On the other hand, when R0=1, one eigenvalue of the expression will be 0. By the concept of Hartman Grobman [24], the ailment-free equilibrium is nonhyperbolic. The absence of non hyperbolicity indicates that the linearized system cannot describe the local features of the system.

4.1 Global Stability of the Ailment-Free Equilibrium

By employing the approach used in [25], the asymptotic stability for the ailment-free equilibrium will be derived in a global sense. This takes the following into account:

ABCDνA′(T)=f(A,y),ABCDνy′(T)=g(A,y),g(A,0)=0,(16)

where A∈ℜm indicates the non-infectious ones and y∈ℜn indicates the infectious ones. Assume U0=(A*,0) represents the ailment-free equilibrium and suppose that:

•   ABCDνA′(T)=f(A,0), A* is globally asymptotically stable,

•   g(A,y)=Xy-g(A,y),g(A,y)≥0 for (A,y)∈X, and X=D1G(A*,0) denotes an M-matrix.

Theorem 4.2 If D1 and D2 are satisfied, and U0 = (A*, 0) is a fixed point, then the model is globally asymptotically stable only if Ro≤1.

Proof. The proof can be shown in a similar way as the process carried out in [24].

5  Parameters Estimation

For the validation of an epidemiological model, it is extremely important to compare the results of simulations with the actual data of infected individuals. This increases the reliability of the proposed disease model. Similar values from the simulations and actual data give better information on the disease being investigated. In addition, unknown values of the working parameters that contribute to the model can be determined.

There are different techniques including maximum likelihood estimation, Bayesian technique, nonlinear least-squares approach, and probability plotting, which can be used to obtain the best parameters. In this research, we utilized the nonlinear least-squares approach for computing the best-fitted parameters, including Π,βh,βe,d,ν,ε, and σ, along with the most important parameter of our proposed model, called the fractional order ν (one of the major components of the study). We also obtained the best-fitted parameters for the classical dysentery system and ABC fractional dysentery system, which are shown in Tab. 1. The actual data for the dysentery diarrhea infected individuals was from Ethiopia, which covered a period of 52 weeks in 2017 [26].

When the least-squares technique is utilized, we need to minimize the objective function. This is achieved by tuning the system’s parameters to fit the available data points accurately. Real data cases for the dysentery diarrhea disease in this study are denoted by m points (xn,yn),n=1,…,m, where xn stands for independent quantity and yn shows dependent quantity. The system function has the structure h(x, r), where s tuned parameters are shown in a vector of parameters, r. Hence, the objective is to identify the parameters which ensure that the system’s simulations for the infected cases fit well with the actual data points.

In the present study, we obtained the best fit by measuring the difference between the real data and the simulations. This is shown below:

En=yn-h(xn,r)(17)

Finally, the optimal set of parameters is obtained, as shown in Tab. 1, with the least-squares approach while minimizing the absolute relative error, on average, as shown below:

ARE=1N∑n=1m|yn-h(xn,r)|yn,(18)

where N stands for the total data value, which is 52 in this study. Moreover, real prevalent cases of the dysentery diarrhea disease, along with the classical and the ABC system’s simulations for the infected individuals, are listed in Tab. 2. In addition, Fig. 1 shows the best fit of the classical and ABC system with the real cases. The ABC system showed an average absolute relative error of 3.4284e −02, and the classical system was 3.4432e −02. Therefore, it shows that the ABC system had some advantages compared to the classical dysentery diarrhea system. In addition, the basic reproductive numbers R0 were 1.7031 and 1.9581 for the ABC and classical systems, respectively. This clearly shows that the disease can be well prevented if the ABC operator is taken into consideration while modeling the epidemic.

Figure 1: Plots and comparison of the real incidence cases, classical (ν=1) and ABC (ν=9.9410e-01) systems for the dysentery diarrhea

6  Sensitivity

In this section, the concept of sensitivity analysis is used to discover the robust significance of the generic parameters present in the base reproduction number R0. Furthermore, both the analytical and numerical values of the R0 parameters are derived from precise assumptions using parameter values.

If the dynamics follow the model (1), the analytical expressions can be used to explain the process of tracking the model’s onset at various locations. The threshold value, R0 can reduce and stop the ailment spread by reducing the number to less than unity. The sensitivity index technique is used to measure the most sensitive parameters in the model. The parameters with a positive sign are considered highly and proportionally sensitive to R0, while those with a negative sign are less sensitive to R0 decreasing. The other category is neutrally sensitive (with zero relative sensitivity). The cause of the transmission of the infringement is directly linked to the specific reproduction number R0. The R0 elasticity indices [27] are shown below:

ΥPiR0=∂R0∂Pi×PiR0,(19)

where R0 denotes the basic reproduction ratio and Pi is stated above. Following the described formula, we get:

Υβh=βhνν(βhνμν+γν+dν+(πβeε)νμνσν(μν+γν+dν)C)-1(μν+γν+dν)-1,Υμ=μ(-βhνμνν(μν+γν+dν)2μ-(πβeε)ννμνσν(μν+γν+dν)Cμ-(πβeε)ννσν(μν+γν+dν)2Cμ)×(βhνμν+γν+dν+(πβeε)νμνσν(μν+γν+dν)C)-1,Υγ=γ(-βhνγνν(μν+γν+dν)2γ-(πβeε)νγννμνσν(μν+γν+dν)2Cγ)×(βhνμν+γν+dν+(πβeε)νμνσν(μν+γν+dν)C)-1,Υd=d(-βhνdνν(μν+γν+dν)2d-(πβeε)νdννμνσν(μν+γν+dν)2Cd)×(βhνμν+γν+dν+(πβeε)νμνσν(μν+γν+dν)C)-1,Υσ=-(πβeε)νν(βhνμν+γν+dν+(πβeε)νμνσν(μν+γν+dν)C)-1(μν)-1(σν)-1(μν+γν+dν)-1C-1,ΥC=-(πβeε)ν(βhνμν+γν+dν+(πβeε)νμνσν(μν+γν+dν)C)-1(μν)-1(σν)-1×(μν+γν+dν)-1C-1,Υπ=(πβeε)νν(βhνμν+γν+dν+(πβeε)νμνσν(μν+γν+dν)C)-1(μν)-1(σν)-1×(μν+γν+dν)-1C-1,Υε=(πβeε)νν(βhνμν+γν+dν+(πβeε)νμνσν(μν+γν+dν)C)-1(μν)-1(σν)-1×(μν+γν+dν)-1C-1,Υβe=(πβeε)νν(βhνμν+γν+dν+(πβeε)νμνσν(μν+γν+dν)C)-1(μν)-1×(σν)-1(μν+γν+dν)-1C-1.(20)

The numerical values indicating the relative significance of R0 are given in Tab. 3. Some parameters are positive while some are negative. Parameters with positive values mean that an increase in the parameter’s values will have a major effect on the frequency of the ailment spread. On the other hand, parameters with negative values mean an increase in such parameters would decrease the effect of the disease. A representation of the values given in Tab. 3, is shown in Fig. 2.

Figure 2: Elasticity indices for various parameters of R0

7  Simulations for the ABC Model

In this section, an algorithm is first developed to obtain the approximate solution of the ABC dysentery diarrhea model, wherein the operator uses non-local and non-singular types of the kernel. The algorithm being developed is discussed in [28], which entails the combination of the fundamental theorem of fractional calculus and two-step Lagrange type polynomial. Therefore, the fundamental theorem of fractional calculus on the Cauchy type initial value problem is given below:

ABCD0νw(t)=G(t,w(t)),w(0)=w0,0<t<∞,(21)

This leads to:

w(t)-w(0)=(1-ν)ABC(ν)G(t,w(t))+νABC(ν)×Γ(ν)∫0tG(τ,w(τ))(t-τ)ν-1dτ.(22)

At t=tn+1,n=0,1,2,…, we have:

w(tn+1)-w(0)=(1-ν)ABC(ν)G(tn,w(tn))+νABC(ν)×Γ(ν)∫0tn+1G(τ,w(τ))(tn+1-τ)ν-1dτ. (23)

w(tn+1)-w(0)=(1-ν)ABC(ν)G(tn,w(tn))+νABC(ν)×Γ(ν)∑v=0n∫tvtv+1G(τ,w(τ))(tn+1-τ)ν-1dτ. (24)

With the help of interpolation polynomial, we approximate function G(τ,w(τ))over[tv,tv+1]:

G(τ,w(τ))≈Pk(τ)=G(tv,w(tv))h(τ-tv-1)-G(tv-1,w(tv-1))h(τ-tv).(25)