# Parameter Estimation and Sensitivity Analysis of Dysentery Diarrhea Epidemic Model

Authors: Hailay Weldegiorgis Berhe; Oluwole Daniel Makinde; David Mwangi Theuri
Journal: Journal of Applied Mathematics
 (2019)

---

## Abstract

In this paper, dysentery diarrhea deterministic compartmental model is proposed. The local and global stability of the disease-free equilibrium is obtained using the stability theory of differential equations. Numerical simulation of the system shows that the backward bifurcation of the endemic equilibrium exists forR0>1. The system is formulated as a standard nonlinear least squares problem to estimate the parameters. The estimated reproduction number, based on the dysentery diarrhea disease data for Ethiopia in 2017, is R0=1.1208. This suggests that elimination of the dysentery disease from Ethiopia is not practical. A graphical method is used to validate the model. Sensitivity analysis is carried out to determine the importance of model parameters in the disease dynamics. It is found out that the reproduction number is the most sensitive to the effective transmission rate of dysentery diarrhea (βh). It is also demonstrated that control of the effective transmission rate is essential to stop the spreading of the disease.


---

## Body

## 1. Introduction

Diarrheal disease is a common threat worldwide, particularly in developing countries. It is preventable and treatable. According to [1], diarrhea in its various forms is usually one of the five major causes of deaths in developing nations. It is the second leading cause of death in children under five years of age. Pneumonia and diarrhea are responsible for an estimated 40 percent of all child deaths each year. The authors [2] also reported that diarrhea causes 18% of all the deaths of children under the age of five. According to the study, more than 5000 children are dying every day as a result of diarrheal diseases. African and South-East Asian regions share 78% of all these deaths.In most cases, diarrhea disease occurs in three forms and frequently occurs in children of under five years age [1, 3]. The first type of diarrhea is called dysentery diarrhea. This is diarrhea with the loose or watery stools with visible blood. Dysentery diarrhea is most often caused by shigella dysenteriae serotype 1 (bacillary dysentery) or Entamoeba histolytica (amoebic dysentery). The second form of diarrhea is acute watery diarrhea. This form includes cholera and is associated with significant fluid loss and rapid dehydration in an infected individual. It usually lasts for several hours or days. This can be caused by the pathogens V. cholerae or E. coli bacteria and rotavirus. The last form of diarrhea is persistent diarrhea. This is diarrhea with or without blood that lasts at least 14 days. It is common in undernourished children and those with other illnesses, such as AIDS.Stochastic and deterministic mathematical models are widely employed by public health practitioners in order to predict and control disease outbreak as well as determining the cost effectiveness of the available strategies [4–10]. In stochastic model, the disease transmission dynamics within a given population is based on the probabilistic analyses, while the deterministic models deals with the use of continuous functions to describe the time evolution of disease in a homogeneously mixed population [11]. Several authors [12–21] have made use of stability theory of differential equations in order to analysis deterministic compartmental models for both local and global dynamics of diseases transmission in a within hosts and between hosts population such as HIV/AIDs, cholera, sheep brucellosis, influenza, typhoid, etc. In all these work, the threshold parameters for disease control and elimination were determined.Epidemiologically speaking, it is important to have adequate knowledge of the physical characteristic of any given ailment or disease in order to provide accurate model that can be utilized to predict and control its outbreak effectively. The availability of relevant data for the disease will also enhance the model parameters estimation and the usefulness of the model with respect to the disease involved. Model parameters from epidemiological models can be easily estimated using methods involving ordinary least squares estimator, maximum likelihood estimator derivative approximation, moments estimator, Markov chain Monte Carlo (MCMC) strategy, derivative-free optimization algorithms, the Levenberg-Marquardt, and Trust-Region-Reflective [22–26]. Biegler and Grossmann [27] employed optimization technique based on the seasonal data with the SIRS epidemic model as a constraint in order to estimate the parameters of a generalized incidence rate function. Li [28] also made use of optimization method with difference equations associated with a given system as a constraint in order to estimate the values of embedded parameters. The SIR model parameters were numerically estimated by Capaldi et al. [29] using least squares method. Other authors [30–34] have made use of various parameter estimation techniques base on the available data to determine the effects of various epidemiological factors on the disease transmission and possible control strategy. However, there is still a need to estimate the parameters involved in the various epidemiological models in order to assess the relative influence of each input variable on the output variable.To the best of our knowledge problem involving estimation of parameters in dysentery diarrhea epidemic model with sensitivity analysis of the basic reproduction number has been scarcely investigated. The main objective of this paper is to tackle this problem by employing constrained optimization technique in order to estimate the values of parameters involved in a deterministic compartmental model for dysentery diarrhea epidemic as a constraint. Real data based on the weekly dysentery diarrhea epidemic that occurred in Ethiopia in 2017 were used to estimate parameters which are not accessible from the characteristic of the disease. In addition, we establish the local and global stability of the disease-free equilibrium and perform the bifurcation analysis. Sensitivity and uncertainty of the basic reproduction number of the system are analyzed. The study estimates the reproduction number of Ethiopia using the estimated parameters to predict the disease spread in the community. In the following Section2, we present the model and its analysis. The numerical approximation and parameter estimation methods are presented in Section 3. In Section 4, numerical simulation of the system and sensitivity analysis with given parameter values is performed. Brief discussion and conclusions are presented in the last section.


## 2. The Model and Its Analysis

The SIRSB (Susceptible-Infected-Recovered-Susceptible-Pathogen population) model for dysentery diarrhea proposed by same authors [35] consists of the following system of nonlinear ordinary differential equations:(1)dSdt=Λ+αR-λh+λB+μS,dIdt=λh+λBS-μ+γ+dI,dRdt=γI-μ+αR,dBdt=ϵI-σB,S0=S0≥0,E0=I0≥0,R0=R0≥0,B0=B0≥0.The standard (frequency dependent) incidence in the human to human interaction and the logistic incidence in the environment to human interaction are respectively given by(2)λh=βhINandλβ=βBBK+B.K is the shigella concentration that yields 50% chance of catching dysentery diarrhea [36]. βB and βh represent rates of ingesting shigella from the contaminated environment and through human to human interaction, respectively. It is assumed that the incidence λB=βBB/K+B is a nonlinear function in B. This incidence represents a Hill function. If B≪K, λB grows linearly with B; if B≫K, λB approaches a constant state, βB, resulting in a saturation state. Therefore, this form is based on the assumption that the incidence has a saturation form. The pathogen population is diminished (at a rate of σ2-σ1=σ>0) [17, 35, 36].The corresponding flow diagram and the description of the parameters for the model are given in Figure1 and Table 1, respectively.Table 1
Description of parameters of the model equation (1).
Parameters Interpretation Units Λ Recruitment rate of susceptible population HumansTime-1 μ Natural death rate of humans Time-1 γ Natural recovery rate of diarrhea Time-1 α Relapse rate of the recovered humans to the susceptible class T i m e - 1 d Disease induced death rate of dysentery diarrhea T i m e - 1 K Concentration of Shigella cells βh Effective transmission rate of diarrhea due to human to human interaction T i m e - 1 βB Effective transmission rate of dysentery diarrhea due to environment to human interaction T i m e - 1 ϵ Pathogen shedding rate CellsHuman-1Time-1 σ1 Shigella Pathogen growth rate T i m e - 1 σ2 Shigella Pathogen death rate T i m e - 1 σ Net death rate of Shigella Pathogen T i m e - 1Figure 1
Schematic diagram for the flow of dysentery diarrhea infection.


### 2.1. Basic Reproduction Number and Existence of Equilibria

System equation (1) has one disease-free equilibrium E0=(Λ/μ,0,0,0). The basic reproduction number (R0) is utilized to quantify the transmission capability of a sickness. It is thought of as the quantity of secondary diseases delivered by a primary case of an infection in a population that is completely susceptible. It can in this manner be estimated by checking the quantity of secondary cases following the presentation of a disease into an absolutely susceptible population. R0 decides if a pathogen can hold on in such a population, and it is important for evaluating control alternatives. At the point when R0 is under 1, on average each infected individual infects less than one other individual, and the pathogen will cease to exist in the population. Conversely, when R0 surpasses 1 there is an exponential ascent in the quantity of cases after some time and disease epidemic results [38]. Using the next generation approach of [39], the reproduction number is(3)R0=βhμ+γ+d+ΛβBϵμμ+γ+dKσ.The endemic equilibrium E∗=(S∗,I∗,R∗,B∗) of the system is the solution of the algebraic equation(4)Λ-λh+λB+μS∗+αR∗=0,λh+λBS∗-μ+d+γI∗=0,γI∗-μ+αR∗=0,ϵI∗-σB∗=0.That is,(5)B∗=ϵI∗σ,R∗=γI∗α+μ,S∗=Λα+μ-μμ+γ+d+αμ+dI∗μα+μandI∗ is the solution of the following quadratic equation:(6)AI∗2+BI∗+C=0,where(7)A=ϵμβh-βBdμμ+d+γ+αμ+d-dμ+d+γ,B=Λϵα+μβBd-μβh+μμ+d+γ+αμ+dKσ+ϵΛβB-μ+d+γKdσ+ϵΛ,C=μΛKσα+μμ+d+γ1-R0.(i)
System equation (1) always has the disease-free equilibrium.(ii)
IfR0>1 and A<0, B>0, then system equation (1) may has two endemic equilibria E1 and E2.(iii)
IfR0>1, and A<0, B<0, then system equation (1) may has a unique endemic equilibrium, E∗.(iv)
IfR0>1 and A>0, B<0, then system equation (1) may has a unique endemic equilibrium, E∗.(v)
IfR0>1 and A>0, B>0, then system equation (1) has no endemic equilibrium.(vi)
If0≤R0<1, there is no endemic equilibrium.


### 2.2. Stability of the Disease-Free Equilibrium

This section is devoted to analytic conditions for the stability of the disease-free equilibrium.


#### 2.2.1. Local Stability of the Disease-Free Equilibrium

Theorem 1.
The disease-free equilibriumE0 is locally asymptotically stable if R0≤1 and is unstable if R0>1.Proof.
Local stability is analyzed by the Jacobian matrix of (1) at E0. That is, (8)JE0=-μ-βhα-βBΛKμ0βh-μ+d+γ0βBΛKμ0γ-α+μ00ϵ0-σ.The eigenvalues of this matrix are(9)λ1=-μ,λ2=-μ+αand  the  solutions  ofλ2+a1λ+a2=0wherea1=μ+d+γ+σ-βhanda2=μ+d+γσ1-R0.If R0<1, a1<0 and a1a2>0. Thus, according to the Hurwitz criterion, the quadratic equation, λ2+a1λ+a2=0, has negative real eigenvalues. Hence, E0 is locally asymptotically stable. If R0=1, one eigenvalue of the quadratic equation is 0 and it is simple. Hence, based on Hartman-Grobman Theorem [40], the disease-free equilibrium E0 is a nonhyperbolic equilibrium. The loss of hyperbolicity of the equilibrium implies that the local behavior of the system cannot be described by the linearized system. And instead center manifold theory has to be used. The bifurcation analysis Theorem 3 below proves that E0 is the only equilibrium at R0=1 and hence E0 is locally asymptotically stable. If R0>1, then E0 is unstable.


#### 2.2.2. Global Stability of the Disease-Free Equilibrium

In this section, we use the method used in [36] to investigate the global asymptotic stability of the disease-free equilibrium. To do so, we write system equation (1) as(10)dXdt=FX,I,dIdt=GX,I,GX,0=0,where X∈Rm denotes the number of noninfectious individuals and I∈Rn denotes the number of infected individuals. Let U0=(X∗,0) denote the disease-free equilibrium of this system. Furthermore, suppose that the following assumptions are true: (H1) For dX/dt=F(X,0), X∗ is globally asymptotically stable;(H2) G(X,I)=AI-G^(X,I), G^(X,I)≥0 for (X,I)∈A,whereA=DIG(X∗,0) is an M-matrix.Theorem 2.
The fixed pointU0=(X∗,0) is a globally asymptotic stable equilibrium of model equation (1) provided that R0≤1 and that assumptions (H1) and (H2) hold.Proof.
To prove this, we write (1) in the form of (10). That is, X=S, I=(I,B). Furthermore, (11)FX,0=Λ-μS0,A=βh-μ+γ+dβBΛKμϵ-σand(12)G^X,I=G1^X,IG2^X,I=1-SNβhI+ΛKμ-SK+BβBB0.Since 1-S/N>0 and Λ/Kμ-S/K+B>0, then G^(X,I)≥0. It is also clear that X∗=(Λ/μ,0) is a global stable equilibrium of the system dX/dt=F(X,0). Hence, the above theorem U0=(X∗,0) is globally stable.


### 2.3. Bifurcation Analysis

Theorem 3.
IfR0<1, the system exhibits a forward bifurcation at R0=1.Proof.
To investigate the type of bifurcation, let(13)βB∗≔μKσϵΛμ+d+γ-βhbe a bifurcation parameter. The Jacobin matrix at E0 with βB=βB∗ is (14)JE0=-μ-βhα-βB∗ΛKμ0βh-μ+d+γ0βB∗ΛKμ0γ-α+μ00ϵ0-σ.The eigenvalues of the characteristics polynomial at R0=1 are λ1=-μ, λ2=-(α+μ), λ3=βhΛ/μ-(μ+d+γ+σ), and λ4=0. It can observe that the three eigenvalues are real and negative and one is 0 and simple. Now we denote by w=(w1,w2,w3,w4)T the right eigenvector corresponding to λ4=0. Direct calculation gives, (15)w=-αμ+d+μμ+d+γμα+μw2,w2,ϵσw2,γμ+αw2T.Moreover, the left eigenvector v=(v1,v2,v3,v4)T that satisfy v.w=1 is given by (16)-μv1=0,-βhv1+βh-μ+d+γv2+γv3+ϵv4=0,αv1-μ+αv3=0,-σϵμ+d+γ-βhv1+σϵμ+d+γ-βhv2-σv4=0.Then, the left eigenvector v turns out to be (17)v=0,v2,0,μ+d+γ-βhϵv2T.The equation v.w=1 results in v2w2+v4w4=1. Assume w2=ϵ(μ+α). Then (18)v=0,1ϵμ+α+γμ+d+γ-βh),0,μ+d+γ-βhϵϵμ+α+γμ+d+γ-βh)T.w=-αμ+d+μμ+d+γμϵ,ϵμ+α,ϵ2σμ+α,γϵT.Based on theoretical results in [40], we have to compute a and b where(19)a=∑k,i,j=13vkwiwj∂2fk∂xi∂xjE0,βB∗,b=∑k,i=13vkwi∂2fk∂xi∂βE0,βB∗.Set x1=S, x2=I, x3=R, x4=B. Evaluating the partial derivatives at (E0,βB∗), we obtain (20)∂2f2∂x1∂x2E0,βB∗=∂2f2∂x2∂x1E0,βB∗=βhS0,∂2f2∂x1∂x4E0,βB∗=∂2f2∂x4∂x1E0,βB∗=-βB∗K∂2f2∂x42E0,βB∗=-2βB∗K2.Furthermore, (21)∂2f2∂x3∂βB∗E0,βB∗=ΛμKand all the other second order partial derivatives are zero. It follows that (22)a=2v2w1w2∂2f2∂x1∂x2E0,βB∗+w1w4∂2f2∂x1∂x4E0,βB∗+w422∂2f2∂x42E0,βB∗and(23)b=v2w4∂2f2∂x4∂βB∗E0,βB∗.Substituting the eigenvectors and the partial derivatives into a and b, we get(24)a=-2ϵμ+αΛϵμ+α+γμ+d+γ-βh)a0whereao=αμ+d+μμ+d+γβh+αμ+d+μ(μ+d+γ+γα+μσϵγα+μμ+d+γ-βhand(25)b=1ϵμ+α+γμ+d+γ-βhγϵΛμK.IfR0<1, then (μ+d+γ)-βh>0, a<0 and b>0. This scenario indicates that the system exhibits transcritical bifurcation at βB=βB∗. Its biological meaning is that as long as R0<1 the dysentery epidemic can be eliminated from the population. If parameters change and result in R0>1, a small endemic state may occur.


## 2.1. Basic Reproduction Number and Existence of Equilibria

System equation (1) has one disease-free equilibrium E0=(Λ/μ,0,0,0). The basic reproduction number (R0) is utilized to quantify the transmission capability of a sickness. It is thought of as the quantity of secondary diseases delivered by a primary case of an infection in a population that is completely susceptible. It can in this manner be estimated by checking the quantity of secondary cases following the presentation of a disease into an absolutely susceptible population. R0 decides if a pathogen can hold on in such a population, and it is important for evaluating control alternatives. At the point when R0 is under 1, on average each infected individual infects less than one other individual, and the pathogen will cease to exist in the population. Conversely, when R0 surpasses 1 there is an exponential ascent in the quantity of cases after some time and disease epidemic results [38]. Using the next generation approach of [39], the reproduction number is(3)R0=βhμ+γ+d+ΛβBϵμμ+γ+dKσ.The endemic equilibrium E∗=(S∗,I∗,R∗,B∗) of the system is the solution of the algebraic equation(4)Λ-λh+λB+μS∗+αR∗=0,λh+λBS∗-μ+d+γI∗=0,γI∗-μ+αR∗=0,ϵI∗-σB∗=0.That is,(5)B∗=ϵI∗σ,R∗=γI∗α+μ,S∗=Λα+μ-μμ+γ+d+αμ+dI∗μα+μandI∗ is the solution of the following quadratic equation:(6)AI∗2+BI∗+C=0,where(7)A=ϵμβh-βBdμμ+d+γ+αμ+d-dμ+d+γ,B=Λϵα+μβBd-μβh+μμ+d+γ+αμ+dKσ+ϵΛβB-μ+d+γKdσ+ϵΛ,C=μΛKσα+μμ+d+γ1-R0.(i)
System equation (1) always has the disease-free equilibrium.(ii)
IfR0>1 and A<0, B>0, then system equation (1) may has two endemic equilibria E1 and E2.(iii)
IfR0>1, and A<0, B<0, then system equation (1) may has a unique endemic equilibrium, E∗.(iv)
IfR0>1 and A>0, B<0, then system equation (1) may has a unique endemic equilibrium, E∗.(v)
IfR0>1 and A>0, B>0, then system equation (1) has no endemic equilibrium.(vi)
If0≤R0<1, there is no endemic equilibrium.


## 2.2. Stability of the Disease-Free Equilibrium

This section is devoted to analytic conditions for the stability of the disease-free equilibrium.


### 2.2.1. Local Stability of the Disease-Free Equilibrium

Theorem 1.
The disease-free equilibriumE0 is locally asymptotically stable if R0≤1 and is unstable if R0>1.Proof.
Local stability is analyzed by the Jacobian matrix of (1) at E0. That is, (8)JE0=-μ-βhα-βBΛKμ0βh-μ+d+γ0βBΛKμ0γ-α+μ00ϵ0-σ.The eigenvalues of this matrix are(9)λ1=-μ,λ2=-μ+αand  the  solutions  ofλ2+a1λ+a2=0wherea1=μ+d+γ+σ-βhanda2=μ+d+γσ1-R0.If R0<1, a1<0 and a1a2>0. Thus, according to the Hurwitz criterion, the quadratic equation, λ2+a1λ+a2=0, has negative real eigenvalues. Hence, E0 is locally asymptotically stable. If R0=1, one eigenvalue of the quadratic equation is 0 and it is simple. Hence, based on Hartman-Grobman Theorem [40], the disease-free equilibrium E0 is a nonhyperbolic equilibrium. The loss of hyperbolicity of the equilibrium implies that the local behavior of the system cannot be described by the linearized system. And instead center manifold theory has to be used. The bifurcation analysis Theorem 3 below proves that E0 is the only equilibrium at R0=1 and hence E0 is locally asymptotically stable. If R0>1, then E0 is unstable.


### 2.2.2. Global Stability of the Disease-Free Equilibrium

In this section, we use the method used in [36] to investigate the global asymptotic stability of the disease-free equilibrium. To do so, we write system equation (1) as(10)dXdt=FX,I,dIdt=GX,I,GX,0=0,where X∈Rm denotes the number of noninfectious individuals and I∈Rn denotes the number of infected individuals. Let U0=(X∗,0) denote the disease-free equilibrium of this system. Furthermore, suppose that the following assumptions are true: (H1) For dX/dt=F(X,0), X∗ is globally asymptotically stable;(H2) G(X,I)=AI-G^(X,I), G^(X,I)≥0 for (X,I)∈A,whereA=DIG(X∗,0) is an M-matrix.Theorem 2.
The fixed pointU0=(X∗,0) is a globally asymptotic stable equilibrium of model equation (1) provided that R0≤1 and that assumptions (H1) and (H2) hold.Proof.
To prove this, we write (1) in the form of (10). That is, X=S, I=(I,B). Furthermore, (11)FX,0=Λ-μS0,A=βh-μ+γ+dβBΛKμϵ-σand(12)G^X,I=G1^X,IG2^X,I=1-SNβhI+ΛKμ-SK+BβBB0.Since 1-S/N>0 and Λ/Kμ-S/K+B>0, then G^(X,I)≥0. It is also clear that X∗=(Λ/μ,0) is a global stable equilibrium of the system dX/dt=F(X,0). Hence, the above theorem U0=(X∗,0) is globally stable.


## 2.2.1. Local Stability of the Disease-Free Equilibrium

Theorem 1.
The disease-free equilibriumE0 is locally asymptotically stable if R0≤1 and is unstable if R0>1.Proof.
Local stability is analyzed by the Jacobian matrix of (1) at E0. That is, (8)JE0=-μ-βhα-βBΛKμ0βh-μ+d+γ0βBΛKμ0γ-α+μ00ϵ0-σ.The eigenvalues of this matrix are(9)λ1=-μ,λ2=-μ+αand  the  solutions  ofλ2+a1λ+a2=0wherea1=μ+d+γ+σ-βhanda2=μ+d+γσ1-R0.If R0<1, a1<0 and a1a2>0. Thus, according to the Hurwitz criterion, the quadratic equation, λ2+a1λ+a2=0, has negative real eigenvalues. Hence, E0 is locally asymptotically stable. If R0=1, one eigenvalue of the quadratic equation is 0 and it is simple. Hence, based on Hartman-Grobman Theorem [40], the disease-free equilibrium E0 is a nonhyperbolic equilibrium. The loss of hyperbolicity of the equilibrium implies that the local behavior of the system cannot be described by the linearized system. And instead center manifold theory has to be used. The bifurcation analysis Theorem 3 below proves that E0 is the only equilibrium at R0=1 and hence E0 is locally asymptotically stable. If R0>1, then E0 is unstable.


## 2.2.2. Global Stability of the Disease-Free Equilibrium

In this section, we use the method used in [36] to investigate the global asymptotic stability of the disease-free equilibrium. To do so, we write system equation (1) as(10)dXdt=FX,I,dIdt=GX,I,GX,0=0,where X∈Rm denotes the number of noninfectious individuals and I∈Rn denotes the number of infected individuals. Let U0=(X∗,0) denote the disease-free equilibrium of this system. Furthermore, suppose that the following assumptions are true: (H1) For dX/dt=F(X,0), X∗ is globally asymptotically stable;(H2) G(X,I)=AI-G^(X,I), G^(X,I)≥0 for (X,I)∈A,whereA=DIG(X∗,0) is an M-matrix.Theorem 2.
The fixed pointU0=(X∗,0) is a globally asymptotic stable equilibrium of model equation (1) provided that R0≤1 and that assumptions (H1) and (H2) hold.Proof.
To prove this, we write (1) in the form of (10). That is, X=S, I=(I,B). Furthermore, (11)FX,0=Λ-μS0,A=βh-μ+γ+dβBΛKμϵ-σand(12)G^X,I=G1^X,IG2^X,I=1-SNβhI+ΛKμ-SK+BβBB0.Since 1-S/N>0 and Λ/Kμ-S/K+B>0, then G^(X,I)≥0. It is also clear that X∗=(Λ/μ,0) is a global stable equilibrium of the system dX/dt=F(X,0). Hence, the above theorem U0=(X∗,0) is globally stable.


## 2.3. Bifurcation Analysis

Theorem 3.
IfR0<1, the system exhibits a forward bifurcation at R0=1.Proof.
To investigate the type of bifurcation, let(13)βB∗≔μKσϵΛμ+d+γ-βhbe a bifurcation parameter. The Jacobin matrix at E0 with βB=βB∗ is (14)JE0=-μ-βhα-βB∗ΛKμ0βh-μ+d+γ0βB∗ΛKμ0γ-α+μ00ϵ0-σ.The eigenvalues of the characteristics polynomial at R0=1 are λ1=-μ, λ2=-(α+μ), λ3=βhΛ/μ-(μ+d+γ+σ), and λ4=0. It can observe that the three eigenvalues are real and negative and one is 0 and simple. Now we denote by w=(w1,w2,w3,w4)T the right eigenvector corresponding to λ4=0. Direct calculation gives, (15)w=-αμ+d+μμ+d+γμα+μw2,w2,ϵσw2,γμ+αw2T.Moreover, the left eigenvector v=(v1,v2,v3,v4)T that satisfy v.w=1 is given by (16)-μv1=0,-βhv1+βh-μ+d+γv2+γv3+ϵv4=0,αv1-μ+αv3=0,-σϵμ+d+γ-βhv1+σϵμ+d+γ-βhv2-σv4=0.Then, the left eigenvector v turns out to be (17)v=0,v2,0,μ+d+γ-βhϵv2T.The equation v.w=1 results in v2w2+v4w4=1. Assume w2=ϵ(μ+α). Then (18)v=0,1ϵμ+α+γμ+d+γ-βh),0,μ+d+γ-βhϵϵμ+α+γμ+d+γ-βh)T.w=-αμ+d+μμ+d+γμϵ,ϵμ+α,ϵ2σμ+α,γϵT.Based on theoretical results in [40], we have to compute a and b where(19)a=∑k,i,j=13vkwiwj∂2fk∂xi∂xjE0,βB∗,b=∑k,i=13vkwi∂2fk∂xi∂βE0,βB∗.Set x1=S, x2=I, x3=R, x4=B. Evaluating the partial derivatives at (E0,βB∗), we obtain (20)∂2f2∂x1∂x2E0,βB∗=∂2f2∂x2∂x1E0,βB∗=βhS0,∂2f2∂x1∂x4E0,βB∗=∂2f2∂x4∂x1E0,βB∗=-βB∗K∂2f2∂x42E0,βB∗=-2βB∗K2.Furthermore, (21)∂2f2∂x3∂βB∗E0,βB∗=ΛμKand all the other second order partial derivatives are zero. It follows that (22)a=2v2w1w2∂2f2∂x1∂x2E0,βB∗+w1w4∂2f2∂x1∂x4E0,βB∗+w422∂2f2∂x42E0,βB∗and(23)b=v2w4∂2f2∂x4∂βB∗E0,βB∗.Substituting the eigenvectors and the partial derivatives into a and b, we get(24)a=-2ϵμ+αΛϵμ+α+γμ+d+γ-βh)a0whereao=αμ+d+μμ+d+γβh+αμ+d+μ(μ+d+γ+γα+μσϵγα+μμ+d+γ-βhand(25)b=1ϵμ+α+γμ+d+γ-βhγϵΛμK.IfR0<1, then (μ+d+γ)-βh>0, a<0 and b>0. This scenario indicates that the system exhibits transcritical bifurcation at βB=βB∗. Its biological meaning is that as long as R0<1 the dysentery epidemic can be eliminated from the population. If parameters change and result in R0>1, a small endemic state may occur.


## 3. Numerical Methods

### 3.1. Parameter Estimation

In this manuscript we solve a dynamics parameter estimation problem. It aims to solve for the unknown parameters of a dynamic model supplied as a system of ordinary differential equations. The problem is setup as a standard nonlinear least squares problem; however the nonlinear function to be fitted involves solving ordinary differential equations using a numerical integration scheme. Least squares are a special case of maximum likelihood in which it is assumed that the data are drawn from a normal distribution with mean given by the solution to the ODE. System equation (1) is formulated as a standard dynamic model in the following form:(26)y′=ft,y,ϑ,yt0=y0,where y is the vector of dependent variables and ϑ is the vector of unknown parameters. The error is sum-of-squares error and is represented by(27)ψϑ=∑i=1yi-yi¯2in which yi¯ is the real data and yi=y(ti,ϑ) is the solution to (26) at time ti for a given ϑ. The aim is to minimize the objective function(28)minϑψϑsubjecttoEq.26to obtain our parameter estimates. The parameter estimation algorithm based on [25] with some modification is summarized below.AlgorithmResult: Optimal estimated parameter values ω.(1)
Guess initial parameter valuesω0. Set ω=ω0.(2)
Using MATLABode45 routine, solve (26) using ω to find the solution yi.(3)
Evaluate error using (27).(4)
Minimize (28) using an optimization algorithmlsqnonlin along withTrust-Region-Reflective to find bounded estimated parameters ω^. Update ω=ω^.(5)
Check convergence criteria. If not converged, go to (2).(6)
On convergence, setω=ω^ to current parameter values.


## 3.1. Parameter Estimation

In this manuscript we solve a dynamics parameter estimation problem. It aims to solve for the unknown parameters of a dynamic model supplied as a system of ordinary differential equations. The problem is setup as a standard nonlinear least squares problem; however the nonlinear function to be fitted involves solving ordinary differential equations using a numerical integration scheme. Least squares are a special case of maximum likelihood in which it is assumed that the data are drawn from a normal distribution with mean given by the solution to the ODE. System equation (1) is formulated as a standard dynamic model in the following form:(26)y′=ft,y,ϑ,yt0=y0,where y is the vector of dependent variables and ϑ is the vector of unknown parameters. The error is sum-of-squares error and is represented by(27)ψϑ=∑i=1yi-yi¯2in which yi¯ is the real data and yi=y(ti,ϑ) is the solution to (26) at time ti for a given ϑ. The aim is to minimize the objective function(28)minϑψϑsubjecttoEq.26to obtain our parameter estimates. The parameter estimation algorithm based on [25] with some modification is summarized below.AlgorithmResult: Optimal estimated parameter values ω.(1)
Guess initial parameter valuesω0. Set ω=ω0.(2)
Using MATLABode45 routine, solve (26) using ω to find the solution yi.(3)
Evaluate error using (27).(4)
Minimize (28) using an optimization algorithmlsqnonlin along withTrust-Region-Reflective to find bounded estimated parameters ω^. Update ω=ω^.(5)
Check convergence criteria. If not converged, go to (2).(6)
On convergence, setω=ω^ to current parameter values.


## 4. Numerical Simulations and Discussion

In this section, we first present the existence and nonexistence of the endemic equilibria. Secondly, fitting the data of dysentery diarrhea cases in Ethiopia is performed to estimate the parameters. Next, the sensitivity of theR0 to the system parameters is analyzed. Finally, the uncertainty analysis is presented using the estimated parameters.


### 4.1. The Disease-Free and Endemic Equilibria

It is proved in Theorem2 that the disease does not exist in the community so long as the reproduction number is less than unity and the disease is endemic in the community if the reproduction number is greater than unity. It is depicted in Figure 2(a). The disease-free equilibrium is globally stable if R0≤1. However, there exists a possibility of multiple equilibria for R0>1 (Figure 2(b)).Figure 2
(a) Forward bifurcation graph forΛ=0.5, βB=0.4, K=1, μ=0.2, ϵ=0.05, d=0.1, γ=0.1, σ=3, α=1, and βh varies βh=0.001:0.01:15. (b) Multiple equilibria for Λ=2, βB=5, K=1, μ=0.5, ϵ=0.05, d=0.8, γ=0.1, σ=3, α=1, and βh varies βh=0:0.001:50.
(a)
Forward bifurcation (b)
Multiple equilibria forR0>1


### 4.2. Dysentery Diarrhea Data

Dysentery diarrhea is one of the most dangerous types of diarrhea which is endemic in Ethiopia. It is reported weekly under Integrated Disease Surveillance and Response System (IDSR) in the Ministry of Health. Ethiopian Weekly Epidemiological Bulletin [41] reports National Surveillance Data Summary of diseases in Ethiopia. For instance, dysentery surveillance data of week 10 of 2018 is 5,345. There is an increment of 614 as compared to week 9, 2018. However, the number of cases reported during this week is lower than the same weeks of the previous two years. The data of dysentery diarrhea reported cases in 2017 used in this simulation are presented in (Table 2). Time series graphs of data given in (Table 2) are shown in Figure 3. As shown in Figure 3, the infected cases have maximum amplitude in the sixteenth week.Table 2
Weekly dysentery diarrhea reports in Ethiopia,2017.
Weeks Infected Weeks Infected Weeks Infected Weeks Infected Weeks Infected 1 4542 12 7000 23 7375 33 7042 43 5583 2 4750 13 7042 24 7167 34 6917 44 5458 3 4792 14 6458 25 7292 35 6208 45 5458 4 5417 15 7083 26 7000 36 5958 46 5458 5 5250 16 7625 27 6583 37 6500 47 5250 6 5125 17 6875 28 6750 38 5958 48 4750 7 5833 18 6375 29 6833 39 5583 49 4750 8 6000 19 6625 30 6792 40 5667 50 4750 9 5542 20 6875 31 6708 41 5667 51 4583 10 5917 21 7000 32 6833 42 5542 52 4500 11 6583 22 7042Figure 3
Weekly time series of dysentery diarrhea data in Ethiopia 2017.


### 4.3. Estimation of Parameters

In this model we have ten parameter values to be estimated. Among these parameters, the value ofμ, natural mortality rate is obtained from the local demography of Ethiopia. The parameter γ, the duration of infectiousness, is also found from dysentery diarrhea characteristics [37]. The initial infectious individuals are taken as the number of infected ones in the first week of 2017. Furthermore, the shigella concentration that yields 25-50% chance of catching dysentery diarrhea is 200 [37]. Thus, the initial parameter values for estimation is ω0=(Λ,βh,βB,K,μ,γ,d,α,ϵ,σ)=(2500,0.0158,0.052,200,0.000457,0.1106,0.0103,0.018, 0.01,0.3319) and the initial population is (S0,I0,R0,B0)=(416253,4542,200,2000). Using the algorithm in Section 3.1 the estimated parameter values are given in Table 3.Table 3
Estimated parameter values.
Parameter Value Bounds Source Λ 462 [ 0,84388.8700 ] Estimated β h 0.69864 [ 0.2175,1.1798 ] Estimated β B 0.04465 [ 0,2.4170 ] Estimated K 200  [36] μ 0.000457  [36] γ 0.1106  [37] d 0.55602 [ 0.0581,1.0540 ] Estimated α 0.0211 [ 0,1.0588 ] Estimated ϵ 0.01  Assumed σ 46.2086 [ 0,2520.3600 ] Estimated S 0 216253 [ 91252,341254 ] Estimated I 0 3152 [ 0,58153 ] Estimated R 0 165 [ 0,25288 ] Estimated B 0 1356 [ 0,4520 ] EstimatedThe reproduction number based on the Ethiopian data is thereforeR0=R2+R1=1.0476+0.0732=1.1208. R2 is the reproduction number due to the human to human interaction and R1 is the reproduction number due to the environment to human interaction.


### 4.4. Sensitivity Analysis

The sensitivity analysis reveals how imperative every parameter is to illness transmission. It is regularly used to decide the robustness of model expectations to parameter values since there are errors in data collection and assumed parameter values [42]. The perturbation of fixed point estimation of model parameters and the uncertainty in the model parameter estimation are the two most commonly used techniques for sensitivity analysis. The sensitivity of a variable with respect to system parameters is usually measured by sensitivity index. Sensitivity indices enable us to quantify the relative change in a variable when a parameter changes. The normalized forward sensitivity index of a variable regarding a parameter is the proportion of the relative change in the variable to the relative change in the parameter. When the variable is a differentiable function of the parameter, the sensitivity index might be on the other hand defined utilizing partial derivatives.Definition 4 (see [42, 43]).
The normalized forward sensitivity index ofR0 that depends differentiably on a parameter ρ is defined by (29)φρR0=∂R0∂ρρR0.For instance,φρR0=1 implies an increase (decrease) of ρ by y% increases (decreases) R0 by the same percentage. On the other hand, φρR0=-1 indicates an increase (decrease) of ρ by y% decreases (increases) R0 by y%. In this case, ρ is called a highly sensitive parameter.Proposition 5.
Given(30)R0=βhμ+γ+d+ΛβBϵμμ+γ+dKσ.Let(31)R2=βhμ+γ+dandR1=ΛβBϵμμ+γ+dKσ.The normalized forward sensitivity index of R0 with respect to the given parameters is given by (32)φβhR0=R2R0,φdR0=-dμ+d+γ,φγR0=-γμ+d+γ,φμR0=-μ2μ+d+γ+R1R0,φβBR0=φϵR0=φΛR0=R1R0,φKR0=φσR0=-R1R0.Proof.
This follows immediately from Definition4.The sensitivity indices ofφΛR0, φβhR0, φβBR0, φϵR0 are positive and the remaining are negative. Since all indices are the functions of other parameters the sensitivity indices will change with the change in other parameter values. The parameter values in Table 3 are used to calculate the sensitivity indices. The normalized sensitivity of the basic reproduction number with respect to the estimated parameters is given in Table 4. For example, φβhR0=0.9901115761; the physical meaning is that increasing (or decreasing) βh by 10% increases (or decreases) R0 by 9.9%. A highly sensitive parameter should be carefully estimated, because a small variation in that parameter will lead to large quantitative changes. As can be seen from Table 4, the basic reproduction number is most sensitive to the effective transmission rate of dysentery diarrhea (βh). The implication of this is that an increase in the transmission rate increases the spread of the disease in the community.Table 4
Estimated parameters.
Parameter Sensitivity indices ofR0 Λ φΛR0=0.009888424079 βh φ β h R 0 = 0.9901115761 βB φ β B R 0 = 0.009888424076 K φ K R 0 = - 0.009888424073 μ φ μ R 0 = - 0.01466434116 γ φ γ R 0 = - 0.1651192337 d φ d R 0 = - 0.8301048495 ϵ φ ϵ R 0 = 0.009888424080 σ φ σ R 0 = - 0.009888424073


### 4.5. Model Validation

Once a model is developed, we use real data to validate it. There are many statistical techniques applied for model validation. The two most frequently used methods are graphical and numerical techniques. Graphical methods demonstrate a comprehensive range of complex features of the relationship between the model and the real data. A numerical technique for model validation, on the other hand, tends to be narrowly targeted on a specific side of the connection between the model and the observed data.


#### 4.5.1. Graphical Method

The residuals from a fitted model are the differences between the real data and the solution to the model. If the model fits the real data appropriately, the residuals approximate the random errors that make the association with the real data and the solution to the system [44]. As a result, if the residuals seem to behave randomly, it is sound to say that the model fits the real data. On the contrary, if nonrandom arrangement is visible in the residuals, the model fails to depict the true representation of the real data. Keeping in mind the goal to check the validity of the model to the given information, Figure 4(a) has been drawn. Figure 4(a) demonstrates the results of fitting the SIRSB model in (1) to the dysentery data in Table 2. For this figure, the continuous line represents the model output and the points represent the real data points. There are a few data points that fit poorly with the model outcome. However, overall the model solution and data agree well.Figure 4
Fitting of the model equation (1) to the reported data in Table 2 and residual analysis using the estimated parameters.
(a)
Fitting (b)
ResidualsResiduals plots for the model and the real data in Table2 are given in Figure 4(b). As given in Figure 4(b), the residuals seem to be random and the standard errors are generally small for the model. In this manner the estimates acquired here are reasonable.The residual analysis is determined by computing autocorrelations for the residuals at varying time Lags. Residual analysis includes two tests: the whiteness test and the independence test. The whiteness test is about autocorrelation function of the residuals for each output. Based on this test, the residual autocorrelation function of a good model is inside the confidence interval of the corresponding estimates, showing that the residuals are uncorrelated. On the other hand, the independence test checks the cross-correlation between the input and the residuals for each input output pair. The independence test makes evident that a good model has the residuals uncorrelated with past inputs. In this case, correlation shows that the model does not describe how part of the output relates to the corresponding input.The autocorrelation function plot for the residuals shows that most of the residuals are not statistically significant (Figure5) which directly follows that the residuals are not highly correlated. Except for one major downturn, the residuals between week one and the last week do not show any particular form. They fluctuate randomly around zero. The independence of the residuals obtained suggests that the parameter estimates obtained are reliable. It is easy to see from the plot that almost all of the autocorrelations except the fifth Lag fall within the 99% confidence limits. One Lag outside the 99% confidence limits does not necessarily show nonrandomness. Furthermore, adjacent observations do not correlate.Figure 5
Correlation of residuals for the infected individualsΛ=462, μ=0.003199, βh=0.69864, βB=0.04465, K=200, α=0.0211, ϵ=0.01, γ=0.1106, σ=46.2086, and K=0.1.


#### 4.5.2. Uncertainties Analysis

In the process of parameter estimation, it is reasonable to assume that the parameters are dependent. Particularly, the parameters could be thought to vary together or covary. Covariance is the measure of how the parameteryi is associated with yj. In the analysis of model parameters a covariance matrix approximates the uncertainties in the model parameters and the possible disorders in the outputs. The diagonal elements of the matrix contain the variances of the parameters and the off-diagonal elements contain the covariances between all potential pairs of parameter. This matrix Σ is positive semidefinite and symmetric. The interrelations between the estimated elements in this are based on the signs. If the covariance between any two estimated parameters is positive, then the parameters values tend to vary in a positive way. If the covariance between any two estimated parameters is 0, then the two parameters are unrelated. On the other hand, if the covariance between any two parameters is negative, then the parameters values tend to move in opposite directions. The covariance matrix Σ for the model is represented below. (33)Σ=βhβBdΛασβh0.230.160.2434654-0.42160βB0.165.630.1552.170.0225869.23d0.240.150.2536864.41-0.45153.64Λ34654521.6536864.407043651412.80-87035.24-762050.23α-0.420.022094-0.45-87035.241.0839.46σ1605869.23153.64-762050.2239.466121430.51.The variance in descending order are recruitment rate of susceptible population, Λ (7043651412.8), the net death rate of shigella pathogen, σ (6121430.51), effective transmission rates of dysentery diarrhea due to the environment to human interaction, βB (5.63), relapse rate of the recovered ones to susceptible, α (1.08), disease induced death rate of dysentery diarrhea, d (0.23), and effective transmission rates of diarrhea due to the human to human interaction, βh (0.23).The covariance between transmission coefficient (βh) and transmission coefficient (βB), recruitment rate of susceptible population (Λ), disease induced death rate of dysentery diarrhea (d), and net death rate of shigella pathogen (σ) is positive. The covariance between transmission coefficient, βB and recruitment rate of susceptible population (Λ), disease induced death rate of dysentery diarrhea (d), net death rate of shigella pathogen (σ), and relapse rate of the recovered ones to susceptible (α) is positive. The covariance between recruitment rate of susceptible population (Λ) and disease induced death rate of dysentery diarrhea (d) is positive. The possible implication of this is that if one parameter is increased, the other parameter must also be increased to achieve a similar fit to the data. On the other hand, the covariance between the relapse rate of the recovered ones to susceptible class (α) and transmission coefficient (βh), disease induced death rate of dysentery diarrhea (d), and recruitment rate of susceptible population (Λ) is negative. The covariance between transmission coefficient (βh) and the relapse rate of the recovered ones to susceptible class (α) is negative. The covariance between recruitment rate of susceptible population (Λ) and the relapse rate of the recovered ones to susceptible class (α) and net death rate of shigella pathogen (σ) is negative. The physical implication of this is that if one parameter is increased, the other parameter must be decreased to achieve a similar fit to the data. Table 5 shows the summary of covariance of the model parameters.Table 5
Covariance of the parameters.
Covariance Direction Covariance Direction Covariance Direction Cβh,βB +ve C β B , d +ve C α , d -ve Cβh,d +ve C β B , Λ +ve C α , Λ -ve Cβh,Λ +ve C β B , α +ve C α , σ +ve Cβh,α -ve C β B , σ +ve C β h , σ +ve C σ , d +ve C σ , Λ -ve


## 4.1. The Disease-Free and Endemic Equilibria

It is proved in Theorem2 that the disease does not exist in the community so long as the reproduction number is less than unity and the disease is endemic in the community if the reproduction number is greater than unity. It is depicted in Figure 2(a). The disease-free equilibrium is globally stable if R0≤1. However, there exists a possibility of multiple equilibria for R0>1 (Figure 2(b)).Figure 2
(a) Forward bifurcation graph forΛ=0.5, βB=0.4, K=1, μ=0.2, ϵ=0.05, d=0.1, γ=0.1, σ=3, α=1, and βh varies βh=0.001:0.01:15. (b) Multiple equilibria for Λ=2, βB=5, K=1, μ=0.5, ϵ=0.05, d=0.8, γ=0.1, σ=3, α=1, and βh varies βh=0:0.001:50.
(a)
Forward bifurcation (b)
Multiple equilibria forR0>1


## 4.2. Dysentery Diarrhea Data

Dysentery diarrhea is one of the most dangerous types of diarrhea which is endemic in Ethiopia. It is reported weekly under Integrated Disease Surveillance and Response System (IDSR) in the Ministry of Health. Ethiopian Weekly Epidemiological Bulletin [41] reports National Surveillance Data Summary of diseases in Ethiopia. For instance, dysentery surveillance data of week 10 of 2018 is 5,345. There is an increment of 614 as compared to week 9, 2018. However, the number of cases reported during this week is lower than the same weeks of the previous two years. The data of dysentery diarrhea reported cases in 2017 used in this simulation are presented in (Table 2). Time series graphs of data given in (Table 2) are shown in Figure 3. As shown in Figure 3, the infected cases have maximum amplitude in the sixteenth week.Table 2
Weekly dysentery diarrhea reports in Ethiopia,2017.
Weeks Infected Weeks Infected Weeks Infected Weeks Infected Weeks Infected 1 4542 12 7000 23 7375 33 7042 43 5583 2 4750 13 7042 24 7167 34 6917 44 5458 3 4792 14 6458 25 7292 35 6208 45 5458 4 5417 15 7083 26 7000 36 5958 46 5458 5 5250 16 7625 27 6583 37 6500 47 5250 6 5125 17 6875 28 6750 38 5958 48 4750 7 5833 18 6375 29 6833 39 5583 49 4750 8 6000 19 6625 30 6792 40 5667 50 4750 9 5542 20 6875 31 6708 41 5667 51 4583 10 5917 21 7000 32 6833 42 5542 52 4500 11 6583 22 7042Figure 3
Weekly time series of dysentery diarrhea data in Ethiopia 2017.


## 4.3. Estimation of Parameters

In this model we have ten parameter values to be estimated. Among these parameters, the value ofμ, natural mortality rate is obtained from the local demography of Ethiopia. The parameter γ, the duration of infectiousness, is also found from dysentery diarrhea characteristics [37]. The initial infectious individuals are taken as the number of infected ones in the first week of 2017. Furthermore, the shigella concentration that yields 25-50% chance of catching dysentery diarrhea is 200 [37]. Thus, the initial parameter values for estimation is ω0=(Λ,βh,βB,K,μ,γ,d,α,ϵ,σ)=(2500,0.0158,0.052,200,0.000457,0.1106,0.0103,0.018, 0.01,0.3319) and the initial population is (S0,I0,R0,B0)=(416253,4542,200,2000). Using the algorithm in Section 3.1 the estimated parameter values are given in Table 3.Table 3
Estimated parameter values.
Parameter Value Bounds Source Λ 462 [ 0,84388.8700 ] Estimated β h 0.69864 [ 0.2175,1.1798 ] Estimated β B 0.04465 [ 0,2.4170 ] Estimated K 200  [36] μ 0.000457  [36] γ 0.1106  [37] d 0.55602 [ 0.0581,1.0540 ] Estimated α 0.0211 [ 0,1.0588 ] Estimated ϵ 0.01  Assumed σ 46.2086 [ 0,2520.3600 ] Estimated S 0 216253 [ 91252,341254 ] Estimated I 0 3152 [ 0,58153 ] Estimated R 0 165 [ 0,25288 ] Estimated B 0 1356 [ 0,4520 ] EstimatedThe reproduction number based on the Ethiopian data is thereforeR0=R2+R1=1.0476+0.0732=1.1208. R2 is the reproduction number due to the human to human interaction and R1 is the reproduction number due to the environment to human interaction.


## 4.4. Sensitivity Analysis

The sensitivity analysis reveals how imperative every parameter is to illness transmission. It is regularly used to decide the robustness of model expectations to parameter values since there are errors in data collection and assumed parameter values [42]. The perturbation of fixed point estimation of model parameters and the uncertainty in the model parameter estimation are the two most commonly used techniques for sensitivity analysis. The sensitivity of a variable with respect to system parameters is usually measured by sensitivity index. Sensitivity indices enable us to quantify the relative change in a variable when a parameter changes. The normalized forward sensitivity index of a variable regarding a parameter is the proportion of the relative change in the variable to the relative change in the parameter. When the variable is a differentiable function of the parameter, the sensitivity index might be on the other hand defined utilizing partial derivatives.Definition 4 (see [42, 43]).
The normalized forward sensitivity index ofR0 that depends differentiably on a parameter ρ is defined by (29)φρR0=∂R0∂ρρR0.For instance,φρR0=1 implies an increase (decrease) of ρ by y% increases (decreases) R0 by the same percentage. On the other hand, φρR0=-1 indicates an increase (decrease) of ρ by y% decreases (increases) R0 by y%. In this case, ρ is called a highly sensitive parameter.Proposition 5.
Given(30)R0=βhμ+γ+d+ΛβBϵμμ+γ+dKσ.Let(31)R2=βhμ+γ+dandR1=ΛβBϵμμ+γ+dKσ.The normalized forward sensitivity index of R0 with respect to the given parameters is given by (32)φβhR0=R2R0,φdR0=-dμ+d+γ,φγR0=-γμ+d+γ,φμR0=-μ2μ+d+γ+R1R0,φβBR0=φϵR0=φΛR0=R1R0,φKR0=φσR0=-R1R0.Proof.
This follows immediately from Definition4.The sensitivity indices ofφΛR0, φβhR0, φβBR0, φϵR0 are positive and the remaining are negative. Since all indices are the functions of other parameters the sensitivity indices will change with the change in other parameter values. The parameter values in Table 3 are used to calculate the sensitivity indices. The normalized sensitivity of the basic reproduction number with respect to the estimated parameters is given in Table 4. For example, φβhR0=0.9901115761; the physical meaning is that increasing (or decreasing) βh by 10% increases (or decreases) R0 by 9.9%. A highly sensitive parameter should be carefully estimated, because a small variation in that parameter will lead to large quantitative changes. As can be seen from Table 4, the basic reproduction number is most sensitive to the effective transmission rate of dysentery diarrhea (βh). The implication of this is that an increase in the transmission rate increases the spread of the disease in the community.Table 4
Estimated parameters.
Parameter Sensitivity indices ofR0 Λ φΛR0=0.009888424079 βh φ β h R 0 = 0.9901115761 βB φ β B R 0 = 0.009888424076 K φ K R 0 = - 0.009888424073 μ φ μ R 0 = - 0.01466434116 γ φ γ R 0 = - 0.1651192337 d φ d R 0 = - 0.8301048495 ϵ φ ϵ R 0 = 0.009888424080 σ φ σ R 0 = - 0.009888424073


## 4.5. Model Validation

Once a model is developed, we use real data to validate it. There are many statistical techniques applied for model validation. The two most frequently used methods are graphical and numerical techniques. Graphical methods demonstrate a comprehensive range of complex features of the relationship between the model and the real data. A numerical technique for model validation, on the other hand, tends to be narrowly targeted on a specific side of the connection between the model and the observed data.


### 4.5.1. Graphical Method

The residuals from a fitted model are the differences between the real data and the solution to the model. If the model fits the real data appropriately, the residuals approximate the random errors that make the association with the real data and the solution to the system [44]. As a result, if the residuals seem to behave randomly, it is sound to say that the model fits the real data. On the contrary, if nonrandom arrangement is visible in the residuals, the model fails to depict the true representation of the real data. Keeping in mind the goal to check the validity of the model to the given information, Figure 4(a) has been drawn. Figure 4(a) demonstrates the results of fitting the SIRSB model in (1) to the dysentery data in Table 2. For this figure, the continuous line represents the model output and the points represent the real data points. There are a few data points that fit poorly with the model outcome. However, overall the model solution and data agree well.Figure 4
Fitting of the model equation (1) to the reported data in Table 2 and residual analysis using the estimated parameters.
(a)
Fitting (b)
ResidualsResiduals plots for the model and the real data in Table2 are given in Figure 4(b). As given in Figure 4(b), the residuals seem to be random and the standard errors are generally small for the model. In this manner the estimates acquired here are reasonable.The residual analysis is determined by computing autocorrelations for the residuals at varying time Lags. Residual analysis includes two tests: the whiteness test and the independence test. The whiteness test is about autocorrelation function of the residuals for each output. Based on this test, the residual autocorrelation function of a good model is inside the confidence interval of the corresponding estimates, showing that the residuals are uncorrelated. On the other hand, the independence test checks the cross-correlation between the input and the residuals for each input output pair. The independence test makes evident that a good model has the residuals uncorrelated with past inputs. In this case, correlation shows that the model does not describe how part of the output relates to the corresponding input.The autocorrelation function plot for the residuals shows that most of the residuals are not statistically significant (Figure5) which directly follows that the residuals are not highly correlated. Except for one major downturn, the residuals between week one and the last week do not show any particular form. They fluctuate randomly around zero. The independence of the residuals obtained suggests that the parameter estimates obtained are reliable. It is easy to see from the plot that almost all of the autocorrelations except the fifth Lag fall within the 99% confidence limits. One Lag outside the 99% confidence limits does not necessarily show nonrandomness. Furthermore, adjacent observations do not correlate.Figure 5
Correlation of residuals for the infected individualsΛ=462, μ=0.003199, βh=0.69864, βB=0.04465, K=200, α=0.0211, ϵ=0.01, γ=0.1106, σ=46.2086, and K=0.1.


### 4.5.2. Uncertainties Analysis

In the process of parameter estimation, it is reasonable to assume that the parameters are dependent. Particularly, the parameters could be thought to vary together or covary. Covariance is the measure of how the parameteryi is associated with yj. In the analysis of model parameters a covariance matrix approximates the uncertainties in the model parameters and the possible disorders in the outputs. The diagonal elements of the matrix contain the variances of the parameters and the off-diagonal elements contain the covariances between all potential pairs of parameter. This matrix Σ is positive semidefinite and symmetric. The interrelations between the estimated elements in this are based on the signs. If the covariance between any two estimated parameters is positive, then the parameters values tend to vary in a positive way. If the covariance between any two estimated parameters is 0, then the two parameters are unrelated. On the other hand, if the covariance between any two parameters is negative, then the parameters values tend to move in opposite directions. The covariance matrix Σ for the model is represented below. (33)Σ=βhβBdΛασβh0.230.160.2434654-0.42160βB0.165.630.1552.170.0225869.23d0.240.150.2536864.41-0.45153.64Λ34654521.6536864.407043651412.80-87035.24-762050.23α-0.420.022094-0.45-87035.241.0839.46σ1605869.23153.64-762050.2239.466121430.51.The variance in descending order are recruitment rate of susceptible population, Λ (7043651412.8), the net death rate of shigella pathogen, σ (6121430.51), effective transmission rates of dysentery diarrhea due to the environment to human interaction, βB (5.63), relapse rate of the recovered ones to susceptible, α (1.08), disease induced death rate of dysentery diarrhea, d (0.23), and effective transmission rates of diarrhea due to the human to human interaction, βh (0.23).The covariance between transmission coefficient (βh) and transmission coefficient (βB), recruitment rate of susceptible population (Λ), disease induced death rate of dysentery diarrhea (d), and net death rate of shigella pathogen (σ) is positive. The covariance between transmission coefficient, βB and recruitment rate of susceptible population (Λ), disease induced death rate of dysentery diarrhea (d), net death rate of shigella pathogen (σ), and relapse rate of the recovered ones to susceptible (α) is positive. The covariance between recruitment rate of susceptible population (Λ) and disease induced death rate of dysentery diarrhea (d) is positive. The possible implication of this is that if one parameter is increased, the other parameter must also be increased to achieve a similar fit to the data. On the other hand, the covariance between the relapse rate of the recovered ones to susceptible class (α) and transmission coefficient (βh), disease induced death rate of dysentery diarrhea (d), and recruitment rate of susceptible population (Λ) is negative. The covariance between transmission coefficient (βh) and the relapse rate of the recovered ones to susceptible class (α) is negative. The covariance between recruitment rate of susceptible population (Λ) and the relapse rate of the recovered ones to susceptible class (α) and net death rate of shigella pathogen (σ) is negative. The physical implication of this is that if one parameter is increased, the other parameter must be decreased to achieve a similar fit to the data. Table 5 shows the summary of covariance of the model parameters.Table 5
Covariance of the parameters.
Covariance Direction Covariance Direction Covariance Direction Cβh,βB +ve C β B , d +ve C α , d -ve Cβh,d +ve C β B , Λ +ve C α , Λ -ve Cβh,Λ +ve C β B , α +ve C α , σ +ve Cβh,α -ve C β B , σ +ve C β h , σ +ve C σ , d +ve C σ , Λ -ve


## 4.5.1. Graphical Method

The residuals from a fitted model are the differences between the real data and the solution to the model. If the model fits the real data appropriately, the residuals approximate the random errors that make the association with the real data and the solution to the system [44]. As a result, if the residuals seem to behave randomly, it is sound to say that the model fits the real data. On the contrary, if nonrandom arrangement is visible in the residuals, the model fails to depict the true representation of the real data. Keeping in mind the goal to check the validity of the model to the given information, Figure 4(a) has been drawn. Figure 4(a) demonstrates the results of fitting the SIRSB model in (1) to the dysentery data in Table 2. For this figure, the continuous line represents the model output and the points represent the real data points. There are a few data points that fit poorly with the model outcome. However, overall the model solution and data agree well.Figure 4
Fitting of the model equation (1) to the reported data in Table 2 and residual analysis using the estimated parameters.
(a)
Fitting (b)
ResidualsResiduals plots for the model and the real data in Table2 are given in Figure 4(b). As given in Figure 4(b), the residuals seem to be random and the standard errors are generally small for the model. In this manner the estimates acquired here are reasonable.The residual analysis is determined by computing autocorrelations for the residuals at varying time Lags. Residual analysis includes two tests: the whiteness test and the independence test. The whiteness test is about autocorrelation function of the residuals for each output. Based on this test, the residual autocorrelation function of a good model is inside the confidence interval of the corresponding estimates, showing that the residuals are uncorrelated. On the other hand, the independence test checks the cross-correlation between the input and the residuals for each input output pair. The independence test makes evident that a good model has the residuals uncorrelated with past inputs. In this case, correlation shows that the model does not describe how part of the output relates to the corresponding input.The autocorrelation function plot for the residuals shows that most of the residuals are not statistically significant (Figure5) which directly follows that the residuals are not highly correlated. Except for one major downturn, the residuals between week one and the last week do not show any particular form. They fluctuate randomly around zero. The independence of the residuals obtained suggests that the parameter estimates obtained are reliable. It is easy to see from the plot that almost all of the autocorrelations except the fifth Lag fall within the 99% confidence limits. One Lag outside the 99% confidence limits does not necessarily show nonrandomness. Furthermore, adjacent observations do not correlate.Figure 5
Correlation of residuals for the infected individualsΛ=462, μ=0.003199, βh=0.69864, βB=0.04465, K=200, α=0.0211, ϵ=0.01, γ=0.1106, σ=46.2086, and K=0.1.


## 4.5.2. Uncertainties Analysis

In the process of parameter estimation, it is reasonable to assume that the parameters are dependent. Particularly, the parameters could be thought to vary together or covary. Covariance is the measure of how the parameteryi is associated with yj. In the analysis of model parameters a covariance matrix approximates the uncertainties in the model parameters and the possible disorders in the outputs. The diagonal elements of the matrix contain the variances of the parameters and the off-diagonal elements contain the covariances between all potential pairs of parameter. This matrix Σ is positive semidefinite and symmetric. The interrelations between the estimated elements in this are based on the signs. If the covariance between any two estimated parameters is positive, then the parameters values tend to vary in a positive way. If the covariance between any two estimated parameters is 0, then the two parameters are unrelated. On the other hand, if the covariance between any two parameters is negative, then the parameters values tend to move in opposite directions. The covariance matrix Σ for the model is represented below. (33)Σ=βhβBdΛασβh0.230.160.2434654-0.42160βB0.165.630.1552.170.0225869.23d0.240.150.2536864.41-0.45153.64Λ34654521.6536864.407043651412.80-87035.24-762050.23α-0.420.022094-0.45-87035.241.0839.46σ1605869.23153.64-762050.2239.466121430.51.The variance in descending order are recruitment rate of susceptible population, Λ (7043651412.8), the net death rate of shigella pathogen, σ (6121430.51), effective transmission rates of dysentery diarrhea due to the environment to human interaction, βB (5.63), relapse rate of the recovered ones to susceptible, α (1.08), disease induced death rate of dysentery diarrhea, d (0.23), and effective transmission rates of diarrhea due to the human to human interaction, βh (0.23).The covariance between transmission coefficient (βh) and transmission coefficient (βB), recruitment rate of susceptible population (Λ), disease induced death rate of dysentery diarrhea (d), and net death rate of shigella pathogen (σ) is positive. The covariance between transmission coefficient, βB and recruitment rate of susceptible population (Λ), disease induced death rate of dysentery diarrhea (d), net death rate of shigella pathogen (σ), and relapse rate of the recovered ones to susceptible (α) is positive. The covariance between recruitment rate of susceptible population (Λ) and disease induced death rate of dysentery diarrhea (d) is positive. The possible implication of this is that if one parameter is increased, the other parameter must also be increased to achieve a similar fit to the data. On the other hand, the covariance between the relapse rate of the recovered ones to susceptible class (α) and transmission coefficient (βh), disease induced death rate of dysentery diarrhea (d), and recruitment rate of susceptible population (Λ) is negative. The covariance between transmission coefficient (βh) and the relapse rate of the recovered ones to susceptible class (α) is negative. The covariance between recruitment rate of susceptible population (Λ) and the relapse rate of the recovered ones to susceptible class (α) and net death rate of shigella pathogen (σ) is negative. The physical implication of this is that if one parameter is increased, the other parameter must be decreased to achieve a similar fit to the data. Table 5 shows the summary of covariance of the model parameters.Table 5
Covariance of the parameters.
Covariance Direction Covariance Direction Covariance Direction Cβh,βB +ve C β B , d +ve C α , d -ve Cβh,d +ve C β B , Λ +ve C α , Λ -ve Cβh,Λ +ve C β B , α +ve C α , σ +ve Cβh,α -ve C β B , σ +ve C β h , σ +ve C σ , d +ve C σ , Λ -ve


## 5. Discussion and Conclusions

In this paper, we have proposed a deterministic compartmental dysentery diarrhea model. We have proved the global stability of the disease-free equilibrium. It is established in Theorem2 that the basic reproduction number R0 is a sharp threshold parameter and it completely determines the global stability of the disease-free equilibrium for R0<1. The possible implication of this is that to eliminate the disease from the community, policymakers have to work on reducing the reproduction number to below unity.The authors [17] have reported that the contribution of each transmission pathway in a disease outbreak is separable to focus on the control efforts of the infectious humans or the contaminated environment. Accordingly, R0 could be written as(34)R0=βhμ+d+γ+ΛβBϵμμ+d+γKσ21-σ1/σ2.Using Taylor’s series expansion, we have(35)1-σ1σ2-1=1+σ1σ2+σ1σ22+higher  order  terms.Thus, (36)R0=R01+R02+R03+higher  order  terms.where R01=βh/μ+d+γ, R02=ϵΛ/μ1/μ+d+γβB/K1/σ2, and R03=ϵΛ/μ1/μ+d+γβB/K1/σ2σ1/σ2.R 01 can be expressed as a primary case in the human population makes an infectious contact with humans at a rate of βh over the mean infectious period of 1/μ+d+γ. R02 can be interpreted as a primary case in the human population makes an infectious contact with the environment at a rate of ϵΛ/μ over the mean infectious period of 1/μ+d+γ and a primary case in the environment makes an infectious contact with humans at a rate of βB/K over the mean infectious period of 1/σ2. R03 is interpreted as a primary case in the human population makes an infectious contact with the environment at a rate of ϵΛ/μ over the mean infectious period of 1/μ+d+γ and a primary case in the environment makes an infectious contact with humans at a rate of βB/K over the mean infectious period of 1/σ2, a fraction σ1/σ2 of which survive and become infectious.Therefore,R0=R01+R02+R03 represents the total contribution to the infectious class and the pathogen population class made by the hosts and pathogens of original case. The higher order terms represent the contribution through direct human to human and environment to human in the higher generation.As an application, we have used our system to simulate the reported dysentery diarrhea cases in Ethiopia in 2017 (Table2) and obtained reasonable matches. To solve the optimization problem, the Runge-Kutta method (ode45) has been used as a solver of a system of nonlinear differential equations along withTrust-Region-Reflective (to find optimum bounded parameter values) andlsqnonlin methods for the sum of squared residuals. The results indicate that simulations of our system can provide a match to the infectious cases in 2017. More importantly, we estimate that the basic reproduction number for dysentery transmission in Ethiopia is R0=1.1208. The implication of this is that dysentery diarrhea is endemic in country. The estimated reproduction number is near to one. If control and treatment measures are implemented, with this value of R0, the disease will be eliminated in short time.The values of parameters describing the system have been estimated by fitting the integrals of the system to the field data on dysentery diarrhea epidemic. Estimation of parameters, in this case, is a challenging task because of missing a large part of the infectious process. One difficulty is that depending on the initially specified parameter values a local minimum may occur and the minimal value of the sum of squared residuals may be different. In such case, the parameter estimation should be repeated with different initial parameter guesses to achieve a better estimate of the global minimum (rather than local) of the sum of squared residuals.



---

*Source: 8465747-2019-02-12.xml*