Mathematical analysis of diarrhoea model with saturated incidence rate

1. Introduction

Diarrhea is the frequent passage of loose, watery, soft stools with or without abdominal bloating, pressure, and cramps commonly referred to as gas or flatulence. It is the second leading cause of death in children under five years old [1]. Diarrhea is responsible for killing around 76,000 children globally, there are nearly 1.7 billion cases of diarrhea disease yearly. In developing countries, the annual incidence rate of diarrhea disease episodes in children less than five years old is 3.2 episodes per child. It kills more young children than HIV/AIDS, malaria and measles combined. It causes more than 1.5 million deaths annually, thereby making it a worse health threat than cancer or AIDS in terms of death toll.

Sub-Sahara / Africa is the most vulnerable region of infectious disease [2], this is due to the fact that the region is greatly affected by climate change which makes it more vulnerable to outbreaks that are associated with periods of rainfall and runoff when subsequent turbidity compromises the efficiency of the drinking water treatment plants [3]. Auld et al. [4] found out that heavy rainfall increases diarrhea outbreak due to water contaminated distribution. Many waterborne disease outbreaks occur following a period of intense rainfall [5]. Diarrhea could be acute which lasts for 2 weeks and chronic which lasts for more than 2 weeks [3]. It is one of the most common diseases that is transferred through contaminated food and water [6,7]. There are two types of diarrhea which are infectious and non-infectious diarrhea.

Infectious diarrhea is caused by a virus, parasite or bacterium, which could be campylobacteria, shiga toxin producing E. Coli, giardiasis, salmonellosis, shigellosis, Rotavirus, yersinia, cryptosphoridiosis etc. Noninfectious is caused by toxins (e.g. food poisoning). This type of diarrhea does not spread from person to person [8, 9]. The immunity after infection is temporary and the infection tends to be less severe than the original infection [8]. However, diarrhea is preventable and can be treated.

Adewale et al. [1] worked on Mathematical analysis of diarrhea in the presence of vaccine. They computed \(R_{0}\). In cases where \(R_{0}>1\), the disease became endemic, meaning the disease remained in the population at a consistent rate, as one infected individual transmits the disease to one susceptible.

Chaturvedi et al. [10] did a study on shigella outbreaks. It was established that as long as \(R_{0}< 1\), there would be no epidemic. Upon simulation using assumed parameter values, the results produced, comprehended the epidemic theory and practical situations. The system was proven stable using the Routh-Hurwitz criterion and parameter estimation was successfully completed.

Jose et al. [11] also worked on Epidemiological model of diarrhea diseases and its application in prevention and control. The model was able to mimic the observed epidemiological patterns of infantile diarrhea diseases associated mainly with enterotoxigenic Escherichia coli or with rotavirus. The proposed mathematical model predicted a plausible pattern of the serological profile of an enteric infection. According to their computer simulation experiments (CSE) with this model, it was not necessary to develop an enteric vaccine conferring total and long-lasting immunity in order to achieve protection from diarrhea diseases in young children. Given a protective efficacy and a finite duration of vaccine-induced protection, the optimal immunization policy must be sought. Oral rehydration therapy (ORT) intervention had a clear effect in diminishing the number of individuals dying from diarrhea illness. The CSE also predict an apparent reduction in age-prevalence of diarrhea diseases by use of ORT.

Ardkaew et al. [12] identified the patterns of diarrhea incidence in children under age five in Northeastern provinces of Thailand along the border with Lao PDR. They based their research on the individual hospital case records of patients with diarrhea from 1999 to 2004. Linear regression models containing the District, seasons, and year as factors were fitted to the log-transformed disease incidences, with generalized estimated equations used to account for spatial correlation between districts. The others observed a higher seasonal trend from January to March and April to June. Their analyses suggested that using a thematic map to show the level of diarrhea incidence by district can help provide practicable information that health authorities can use to work effectively and initiate health policies to eradicate the disease.

Cherry et al. [5] worked on Evaluation of bovine viral diarrhea virus control using a mathematical model of infection dynamics. The model architecture was a development of the traditional model framework using susceptible, infectious and removed animals (the SIR model). The model predicted 1:2% persistent infection (within the range of field estimates) and was fairly insensitive to alterations of structure or parameter values. This model drew important conclusions regarding the control of BVD, particularly with respect to the importance of persistently infected (PI) animals in maintaining BVD as an endemic entity in the herd.

Diarrhea modeling has increasingly attracted the attention of mathematical modelers and several works in that regards have been considered on this disease [13]. However, most of these works have focused on vaccination. Thus, often vaccination of diarrhea is considered for children and adults in susceptible class and the exposed periods. Being motivated by the work done in [14], we included saturation treatment function to get its effect. The aim of this work is to construct a mathematical model to study the transmission dynamics of diarrhea infections within the population.

2. Model Formulation

In this section, we formulate a compartmental model of diarrhea with saturated incidence rate to gain insight of the effect of saturation function on treatment of diarrhea disease.

2.1. Model Diagram

This model subdivides into four (4) compartments namely susceptible class (S), infective class (I), treated class (T) and recovered class (R). Thus, at time \(t\), the total population: \begin{equation}\tag{A} N(t)=S(t)+I(t)+T(t)+R(t) \end{equation} In Table 1 we give the detailed explanation of the parameters and variables used in model A

Table 1. Description of Variables and Parameters of the model.

Parameters	Description
\(S\)	Susceptible
\(I\)	Infective
\(T\)	Treated class
\(R\)	Recovered class
\(\Lambda \)	Recruitment rate
\(\beta _{1} \)	Effective contact rate
\(\beta _{2} \)	Saturation treatment rate
\(\eta \)	Enhancement factor
\(P\)	Proportion of infected individuals joining either the class R or T
\(\gamma \)	Rate of treated individuals from infection
\(\sigma \)	Rate at which treated individuals move to recovered class
\(\omega \)	Rate at which recovered individuals move to susceptible class
\(\mu \)	The natural death rate in all compartments
\(\kappa \)	Educational adjustment

2.2. Model Equations

The following nonlinear system of dynamical equations was derived from model A: Subjected to non-negative initial conditions \(S=S_{0} \ge 0,I=I_{0} \ge 0,T=T_{0} \ge =,R=R_{0} \).

2.3. Model analysis

In this section, we find the fundamental properties of the system (1), which is essential in the proof of the proceeding sections.

2.3.1. Positivity and Boundedness

The associated parameters of the system (1) with respect to the initial conditions are non-negative for all \(t>0\) and we prove this in the following procedures to get the results.

Lemma 1. If \(\left\{S(0),I(0),T(0),R(0)\right\}\) and all associated parameters of the system are positive, then solutions \(S(t),I(t),T(t),R(t)\) are all positive \(\forall t>0\).

Proof. Let \[t_{1} =\sup \left\{t>0:S(t)\ge 0,I(t)\ge 0,T(t)\ge 0,R(t)\ge 0\right\}\] From the first equation of (1), we have \[\frac{dS}{dt} =\Lambda +\omega R-\mu S-\left(\beta _{1} \frac{I}{1+\kappa I} +\beta _{2} \frac{\eta T}{1+\kappa T} \right)S.\] Since \(\left(\Lambda +\omega R\right)\ge 0\), then \(\frac{dS}{dt} \ge -\left(\mu +\phi \right)S,\)
where \(\phi =\left(\beta _{1} \frac{I}{1+\kappa I} +\beta _{2} \frac{\eta T}{1+\kappa T} \right)\).
Integrating both sides gives \[\int _{0}^{t_{1} }\frac{dS}{dt} \ge -\int _{0}^{t_{1} }\left(\mu +\phi \right)dt ,\] \[\left. \left(\ln \left|S(t)\right|\right)\right|_{0}^{t_{1} } \ge -\left(\mu +\phi \right)t_{1} ,\] \[\left. S(t)\right|_{0}^{t_{1} } \ge e^{-\left(\mu +\phi \right)t_{1} } ,\] \[S(t)\ge S(0)+e^{-\left(\mu +\phi \right)t_{1} } >0.\] Next, we want to prove that \(I(t)>0\). From the second equation of (1), we have

Since \(\phi S\ge 0\), then we get \[\frac{dI}{dt} \ge -\left(\mu +\gamma \right)I\] Using the approach from (2), we get \[\int _{0}^{t_{1} }\frac{dI}{I} \ge -\int _{0}^{t_{1} }\left(\mu +\gamma \right)dt ,\] \[\left. \left(\ln \left|I(t)\right|\right)\right|_{0}^{t_{1} } \ge -\left(\mu +\gamma \right)t_{1} ,\] \[\left. I(t)\right|_{0}^{t_{1} } \ge e^{-\left(\mu +\gamma \right)t_{1} } ,\] \[I(t_{1} )\ge I(0)+e^{-\left(\mu +\gamma \right)t_{1} } >0.\] Also we prove that \(T(t)>0\). From the third equation of (1), we have \[\frac{dT}{dt} =P\gamma I-\left(\mu +\sigma \right)T\] Since \(P\gamma I\ge 0\) then we have \[\frac{dT}{dt} \ge -\left(\mu +\sigma \right)T. \] This yields \[\int _{0}^{t_{1} }\frac{dT}{T} \ge -\int _{0}^{t_{1} }\left(\mu +\sigma \right)dt \] \[\left. \left(\ln \left|T(t)\right|\right)\right|_{0}^{t_{1} } \ge -\left(\mu +\sigma \right)t_{1} \] \[\left. T(t)\right|_{0}^{t_{1} } \ge e^{-\left(\mu +\sigma \right)t_{1} } ,\] \[T(t_{1} )\ge T(0)+e^{-\left(\mu +\sigma \right)t_{1} } >0.\] Again, we prove that \(R(t)\ge 0\). From the fourth equation of Equation (1) we have \[\frac{dR}{dt} =\left(1-P\right)\gamma I+\sigma T-\left(\mu +\omega \right)R.\] Since \(\left(1-P\right)\gamma I+\sigma T\ge 0\). Then, we have \[\frac{dR}{dt} \ge -\left(\mu +\omega \right)R. \] Integrating gives us \[\int _{0}^{t_{1} }\frac{dR}{R} \ge -\int _{0}^{t_{1} }\left(\mu +\omega \right)dt ,\] \[\left. \left(\ln \left|R(t)\right|\right)\right|_{0}^{t_{1} } \ge -\left(\mu +\omega \right)t_{1} ,\] \[\left. R(t)\right|_{0}^{t_{1} } \ge e^{-\left(\mu +\omega \right)t_{1} } ,\] \[R(t)\ge R_{0} +e^{-\left(\mu +\omega \right)t_{1} } >0.\] Also, the rate at which the total population varies over time is given by:

Lemma 2. The closed set \(\Omega =\left\{\left(S,I,T,R\right)\in {R}_{+}^{4} :0\le \left(S+I+T+R\right)\le \frac{\Lambda }{\mu } \right\}\) is positively invariant.

Proof. Assume that \(\left\{S(t),I(t),T(t).R(t)\right\}\in {R}_{+}^{4} \forall t>0\). The equation (3) can be written as \[\frac{dN}{\Lambda -\mu N} =dt\] \[\int _{o}^{t}\frac{dN}{\Lambda -\mu N} =\int _{0}^{t}dt \] \[N(t)=N(0)e^{-\mu t} +\frac{\Lambda }{\mu } \left(1-e^{-\mu t} \right)\] \[{\mathop{\lim }\limits_{t\to \infty }} N(t)=\frac{\Lambda }{\mu } \] if \(N(0)\le \frac{\Lambda }{\mu }\), then we have \(N(t)\le \frac{\Lambda }{\mu } ,\forall t>0\). Moreover, if \(N(0)> \frac{\Lambda }{\mu }\), then the solution \(\left(S(t),I(t),T(t),R(t)\right)\) enter the closed set \(\Omega \) which affirms that \(\Omega \) is positively invariant. So, the region \(\Omega\) contains all solutions in \({R}_{+}^{4}\). Hence, it is sufficient to study the disease transmission dynamics under the dynamical system (1) in \(\Omega\).

3. Stability Analysis

In this section, we determine the disease free and endemic equilibrium points of Equation (1).

3.1. Disease Free Equilibrium (DFE)

The disease-free state denoted \(E^{0}\) is when there is no infection i.e \(I=T=R=0\) and is obtained by taking the right side of Equation (1) equal to zero. The corresponding disease-free equilibrium point is \[E^{0} =\left(S^{0} ,0,0,0\right)=\left(\frac{\Lambda }{\mu } ,0,0,0\right).\]

3.2. Endemic Equilibrium (EE)

In this subsection, we determine the EE of Equation (1) denoted \(E^{1} =\left(S^{*} ,I^{*} ,T^{*} ,R^{*} \right)\). Taking the right side of Equation (1) equal to zero, the corresponding EE point is given by the following process. Solving the equations, \[T^{*} =\frac{P\gamma i_{h} }{\mu +\sigma } ,\] \[R^{*} \frac{\gamma \mu i_{h} -P\gamma \mu i_{h} +\gamma \sigma i_{h} }{(\mu +\sigma )\left(\mu +\omega \right)} ,\] \[S^{*} =\frac{-\Lambda +\frac{\omega \left(\gamma \mu i_{h} -P\gamma \mu i_{h} +\gamma \sigma i_{h} \right)}{(\mu +\sigma )\left(\mu +\omega \right)} }{-\mu -\frac{i_{h} \beta _{1} }{1+ki_{h} } -\frac{P\gamma \eta i_{h} \beta _{2} }{(\mu +\sigma )\left(1+\frac{kP\gamma i_{h} }{\mu +\sigma } \right)} } =\odot (i_{h} ).\] Finally by substituting \(S^{*} ,T^{*} \) in the equation below \[\frac{I_{h} \beta _{1} S^{*} }{\left(1+kI_{h} \right)} +\frac{\beta _{2} \eta T^{*} S^{*} }{\left(1+kT^{*} \right)} -\left(\mu +\gamma \right)I_{h} =0,\] Where \(i_{h}\) stands for \(I_{h}^{*}\) and simplifying; the following polynomial equation of order three holds. That is \(i_{h} \left(a_{0} i_{h}^{2} +a_{1} i_{h} +a_{2} \right)=0\) and the coefficients’ \(\left(a_{0} ,a_{1} ,a_{2} \right)\) expressions are given below \begin{eqnarray*} a_{0}&=&kP\gamma \left(k\mu \left(\gamma +\mu \right)\left(\mu +\sigma \right)\left(\mu +\omega \right)+\left(\mu \right)\left(\mu +\sigma \right)\left(\mu +\omega \right)+\gamma \left(\mu ^{2} +2\sigma \omega +\mu \left(\sigma +2\omega -P\omega \right)\right)\right)\left(\beta _{1} +\eta \beta _{2} \right),\\ a_{1} &=&-k\mu \left(\gamma +\mu \right)\left(\mu +\sigma \right)\left(P\gamma +\mu +\sigma \right)\left(\mu +\omega \right)-\left(\mu +\sigma \right)\left(\mu \left(\gamma +\mu \right)\left(\mu +\sigma \right)\right)\\ &&+\left(\mu \left(\mu +\sigma \right)+\gamma \left(2\mu -P\mu +2\sigma \right)\right)\omega -kP\gamma \Lambda \left(\mu +\omega \right)\beta _{1} \\ &&+P\gamma \eta -\gamma \mu \left(\mu +\sigma \right)+\gamma \left(\left(-2+P\right)\mu -2\sigma \right)\omega \\ &&+\left(k\Lambda -\mu \right)\left(\mu +\sigma \right)\left(\mu +\omega \right)\beta _{2}\\ a_{2} &=&-\left(\mu +\sigma \right)\left(\mu +\omega \right)[1-R_{0}] \end{eqnarray*} Therefore, the endemic equilibrium point is given by \[E^{1} =\left(S^{*} ,I^{*} ,T^{*} ,R^{*} \right)\] \[=\left(\odot \left(i_{h} \right),i_{h} ,\frac{P\gamma i_{h} }{\mu +\sigma } ,\frac{\gamma \mu i_{h} -P\gamma \mu i_{h} +\gamma \sigma i_{h} }{\left(\mu +\sigma \right)\left(\mu +\omega \right)} \right).\]

3.3. Basic Reproduction Number

The basic reproduction number denoted by \(R_{0}\) is defined as the average number of secondary infection generated by one infected individual into a completely susceptible population at a time t. It is a threshold parameter that allows us to predict whether the disease will die out or persist [11]. We calculate our \(R_{0}\) using the next generation matrix [15]. By referring to this approach, it is given by \(R_{0} =\rho \left(FV^{-1} \right)\) where \(\rho \left(A\right)\) denotes the spectral radius of a matrix \(A, F\) is new infection transfer terms, \(V\) is the non-singular matrix of the remaining transfer terms. \[F=\left[\begin{array}{c} {\beta _{1} \frac{IS}{1+kI} +\beta _{2} \frac{\eta TS}{1+kT} } \\ {0} \end{array}\right].\] Then the Jacobian of \(F\) at \(E^{0} \) is \[F=\left[\begin{array}{cc} {\frac{\Lambda \beta _{1} }{\mu } } & {\frac{\eta \Lambda \beta _{2} }{\mu } } \\ {0} & {0} \end{array}\right].\] Let \(V\) denotes the Jacobian of \(v\) at \(E^{0} \), \[v=\left[\begin{array}{c} {\left(\mu +\gamma \right)I} \\ {-P\gamma I+\left(\mu +\sigma \right)T} \end{array}\right].\] \[V=\left[\begin{array}{cc} {\left(\mu +\gamma \right)} & {0} \\ {-P\gamma } & {\left(\mu +\sigma \right)} \end{array}\right].\] So, below is the next generation matrix obtained, and it allows us to find the basic reproduction number \[A=FV^{-1} =\left[\begin{array}{cc} {\frac{\Lambda \beta _{1} }{\mu \left(\mu +\gamma \right)} +\frac{p\gamma \eta \Lambda \beta _{2} }{\mu \left(\mu +\gamma \right)\left(\mu +\sigma \right)} } & {\frac{\eta \Lambda \beta _{2} }{\mu \left(\mu +\sigma \right)} } \\ {0} & {0} \end{array}\right].\] \(R_{0} \) is the leading eigenvalue of the next generation matrix, \[R_{0} =\frac{\Lambda \left(\mu \beta _{1} +\sigma \beta _{1} +P\gamma \eta \beta _{2} \right)}{\mu \left(\mu +\gamma \right)\left(\mu +\sigma \right)} .\]

3.4. Local Stability of Disease-Free Equilibrium (\(E^{0} \))

In this subsection, we prove the local stability of the disease-free equilibrium of Equation (1) .

Theorem 3. The disease-free equilibrium is locally asymptotically stable if \(R_{0} < 1\). \[J=\left[\begin{array}{cccc} {-\left(\mu +\phi \right)} & {-\beta _{1} \frac{S\left(1+kI\right)-ISk}{\left(1+kI\right)^{2} } } & {-\beta _{2} \frac{\eta S\left(1+kT\right)-\eta TSk}{\left(1+kT\right)^{2} } } & {\omega } \\ {\phi } & {-\beta _{1} \frac{S\left(1+kI\right)-ISk}{\left(1+kI\right)^{2} } -\left(\mu +\gamma \right)} & {-\beta _{2} \frac{\eta S\left(1+kT\right)-\eta TSk}{\left(1+kT\right)^{2} } } & {0} \\ {0} & {P\gamma } & {-\left(\mu +\sigma \right)} & {0} \\ {0} & {\left(1-P\right)\gamma } & {\sigma } & {-\left(\mu +\omega \right)} \end{array}\right],\] where \(\phi =\left(-\beta _{1} \frac{I}{\left(1+kI\right)} +-\beta _{1} \frac{\eta T}{\left(1+kI\right)} \right)\).

Evaluating the Jacobian at \(E_{0}\). Hence, we obtain \[J_{E^{0} } =\left[\begin{array}{cccc} {-\mu } & {-\beta _{1} \frac{\Lambda }{\mu } } & {-\beta _{2} \frac{\eta \Lambda }{\mu } } & {\omega } \\ {\phi } & {\beta _{1} \frac{\Lambda }{\mu } -\left(\mu +\gamma \right)} & {\beta _{2} \frac{\eta \Lambda }{\mu } } & {0} \\ {0} & {P\gamma } & {-\left(\mu +\sigma \right)} & {0} \\ {0} & {\left(1-P\right)\gamma } & {\sigma } & {-\left(\mu +\omega \right)} \end{array}\right].\] Here, we evaluate the Jacobian at \(E^{0}\) and identity matrix to get the eigenvalues, that is, \[\left|J_{E^{0} } -\lambda I\right|=\left|\begin{array}{cccc} {-\mu -\lambda } & {-\beta _{1} \frac{\Lambda }{\mu } } & {-\beta _{2} \frac{\eta \Lambda }{\mu } } & {\omega } \\ {\phi } & {\beta _{1} \frac{\Lambda }{\mu } -\left(\mu +\gamma \right)-\lambda } & {\beta _{2} \frac{\eta \Lambda }{\mu } } & {0} \\ {0} & {P\gamma } & {-\left(\mu +\sigma \right)-\lambda } & {0} \\ {0} & {\left(1-P\right)\gamma } & {\sigma } & {-\left(\mu +\omega \right)-\lambda } \end{array}\right|=0.\] Then \(\lambda _{1} =-\mu < 0\) and \(\lambda _{2} =-\left(\mu +\omega \right)< 0\).The remaining eigenvalues can be obtained by solving \[H=\left|\begin{array}{cc} {\beta _{1} \frac{\Lambda }{\mu } -\left(\mu +\gamma \right)-\lambda } & {\beta _{2} \frac{\eta \Lambda }{\mu } } \\ {P\gamma } & {-\left(\mu +\sigma \right)-\lambda } \end{array}\right|=0.\] Then we get the following characteristic equation \[\lambda ^{2} +\left(2\mu +\gamma +\sigma -\beta _{1} \frac{\Lambda }{\mu } \right)\lambda +\mu ^{2} +\left(\sigma +\gamma -\beta _{1} \frac{\Lambda }{\mu } \right)\mu +\left(\gamma -\beta _{1} \frac{\Lambda }{\mu } \right)-\beta _{2} \frac{P\gamma \eta \Lambda }{\mu } =0.\] We have a form of quadratic equation, that is, \[P\left(\lambda \right)=\lambda ^{2} +a_{1} \lambda +a_{2} =o,\] where \(a_{1} =\left(\mu +\gamma \right)+\left(\mu +\sigma \right)-\beta _{1} \frac{\Lambda }{\mu }\) and \(a_{2} =\left(\mu +\sigma \right)\left(\mu +\gamma \right)\left[1-R_{0} \right]\). If \(R_{0} 0\) and if \(R_{0} < 1\), then \(\beta _{1} \frac{\Lambda }{\mu \left(\mu +\gamma \right)} 0\). Hence the disease-free equilibrium is locally asymptotically stable if \(R_{0} < 1\).

3.5. Global Stability of the Disease Free Equilibrium (\(E^{0}\))

Theorem 4. The disease free equilibrium is globally stable in \(\Omega\) if \(R_{0} < 1\), where \(\Omega\) is a feasible region of Equation (1) and contains all solutions in \({R}_{+}^{4}\) as shown in the previous proof. Again, we prove that the disease-free equilibrium is globally asymptotically stable, using the approach of Castillo – Chavez [13].

Proof. We group Equation (1) into two compartments, that is, uninfected and infected individuals, given by \[\left. \begin{array}{l} {a_{1} :\frac{dX}{dt} =F(X,Z)} \\ {a_{2} :\frac{dZ}{dt} =G(X,Z),G(X,0)=0} \end{array}\right\}. \] Where \(X=(S,R),Z=(I,T)\) with \(X\in {R}_{+}^{2}\) representing the uninfected individuals and \(Z\in {R}_{+}^{2}\) representing infected individuals including the latent and infectious. Let us denote the disease-free equilibrium point by

The point \(E^{0} =(N_{0} ,0)\) is globally asymptotically stable equilibrium of Equation (5) if the following conditions are satisfied: \(b_{1} :E^{0}\) is globally asymptotically stable for \(\frac{dX}{dt} =F(X,0)\), \[b_{2} :\hat{G}(X,Z)\ge 0,(X,Z)\in \Omega ,\] Where \(G(X,Z)=AZ-\hat{G}(X,Z),A=D_{Z} G(N_{0} ,0)\) is a Metzler matrix. Then we can write the matrix \(A\)in the form of \(A=F-V\) and \(\Omega\) is given by Equation (5). Then we have: \[\frac{dX}{dt} =F(X,Z)=\left[\begin{array}{c} {\Lambda +\omega R-\left(\phi +\mu \right)S} \\ {\left(1-P\right)\gamma I+\sigma T-\left(\mu +\omega \right)R} \end{array}\right].\] where \(\phi =\left(\beta _{1} \frac{I}{1+kI} +\beta _{2} \frac{\eta T}{1+kT} \right).\) \[\Rightarrow F(X,0)=\left[\begin{array}{c} {\Lambda -\mu S} \\ {0} \end{array}\right].\] \(\frac{dZ}{dt} =G(X,Z)=\left[\begin{array}{c} {\phi S-\left(\mu +\gamma \right)I} \\ {P\gamma I-\left(\mu +\sigma \right)T} \end{array}\right]\), and \(G(X,0)=0\). Hence \(N_{0} =\left(\frac{\Lambda }{\mu } ,0\right)\) is globally asymptotically stable equilibrium of \(\frac{dX}{dt} =F(X,0)\), that is, \[b_{1} :\frac{dX}{dt} =F(X,0)=\left[\begin{array}{c} {\Lambda -\mu S} \\ {0} \end{array}\right].\] Solving for \(b_{2}\), yields \[A=F-V=\left[\begin{array}{cc} {\beta _{1} \frac{\Lambda }{\mu } -\left(\mu +\gamma \right)} & {\beta _{2} \frac{\eta \Lambda }{\mu } } \\ {P\gamma } & {-\left(\mu +\sigma \right)} \end{array}\right],\] with \(AZ\) given by: \[AZ=\left[\begin{array}{cc} {\beta _{1} \frac{\Lambda }{\mu } -\left(\mu +\gamma \right)} & {\beta _{2} \frac{\eta \Lambda }{\mu } } \\ {P\gamma } & {-\left(\mu +\sigma \right)} \end{array}\right]\left[\begin{array}{c} {I} \\ {T} \end{array}\right]=\left[\begin{array}{c} {\beta _{1} \frac{\Lambda }{\mu } I-\left(\mu +\gamma \right)I+\beta _{2} \frac{\eta \Lambda T}{\mu } } \\ {P\gamma I-\left(\mu +\sigma \right)T} \end{array}\right].\] Thus, \[\hat{G}\left(X,Z\right)=AZ-G\left(X,Z\right)=\left[\begin{array}{c} {0} \\ {0} \end{array}\right]=\left[\begin{array}{c} {\hat{G}\left(X,Z\right)_{1} } \\ {\hat{G}\left(X,Z\right)_{2} } \end{array}\right].\] Since \(\hat{G}\left(X,Z\right)_{1} =0,\) and \(\hat{G}\left(X,Z\right)_{2} =0\), we have \(\hat{G}\left(X,Z\right)\ge 0\) for \(\left(X,Z\right)\in \Omega\). Therefore, the disease-free equilibrium \(E^{0}\) is globally asymptotically stable if \(R_{0} <1 \).

3.6. Local Stability of Endemic Equilibrium

In this subsection we discuss about the local stability of the endemic equilibrium state of Equation (1) .

Theorem 5. The endemic equilibrium \(E^{1}\) is locally asymptotically stable if \(R_{0} >1\). Now, evaluating the Jacobian at \(E^{1}\), we obtain \[J_{E^{1} } =\left[\begin{array}{cccc} {-\left(\mu +A\right)} & {-B} & {-C} & {\omega } \\ {A} & {B-D} & {C} & {0} \\ {0} & {P\gamma } & {-F} & {0} \\ {0} & {\left(1-P\right)\gamma } & {\sigma } & {-U} \end{array}\right],\] Where \(A=\left(\beta _{1} \frac{I^{*} }{1+kI^{*} } +\beta _{2} \frac{\eta T^{*} }{1+kT^{*} } \right),\) \[B=\beta _{1} \frac{S^{*} \left(1+kI^{*} \right)-IS^{*} k}{\left(1+kI^{*} \right)^{2} } ,\] \[C=\beta _{2} \frac{\eta S^{*} \left(1+kT^{*} \right)-\eta T^{*} S^{*} k}{\left(1+kT^{*} \right)^{2} } ,\] \[D=\left(\mu +\gamma \right),\] \[F=\left(\mu +\sigma \right),\] \[U=\left(\mu +\omega \right).\] The eigenvalues of the Jacobian matrix above are, \[\left|J_{E^{1} } -\lambda I\right|=\left|\begin{array}{cccc} {-\left(\mu +A\right)-\lambda } & {-B} & {-C} & {\omega } \\ {A} & {B-D-\lambda } & {C} & {0} \\ {0} & {P\gamma } & {-F-\lambda } & {0} \\ {0} & {\left(1-P\right)\gamma } & {\sigma } & {-U-\lambda } \end{array}\right|=0\] This has the following form of characteristic equation \[-\left(\left(\mu +A\right)+\lambda \right)\left[-\left(U+\lambda \right)\right]\left(\lambda ^{2} +a_{1} \lambda +a_{2} \right)=0.\] Where \begin{eqnarray*} a_{1} &=&-k\mu \left(\gamma +\mu \right)\left(\mu +\sigma \right)\left(P\gamma +\mu +\sigma \right)\left(\mu +\omega \right)-\left(\mu +\sigma \right)\left(\mu \left(\gamma +\mu \right)\left(\mu +\sigma \right)\right)\\ &&+\left(\mu \left(\mu +\sigma \right)+\gamma \left(2\mu -P\mu +2\sigma \right)\right)\omega -kP\gamma \Lambda \left(\mu +\omega \right)\beta _{1}\\ &&+P\gamma \eta \left(-\gamma \mu \left(\mu +\sigma \right)+\gamma \left(\left(-2+P\right)\mu -2\sigma \right)\omega +\left(k\Lambda -\mu \right)\left(\mu +\sigma \right)\left(\mu +\omega \right)\right)\beta _{2} \end{eqnarray*} and \[a_{2} =-\left(\mu +\sigma \right)\left(\mu +\omega \right)\left[1-R_{0} \right].\]

We can see that \(\lambda _{1} =-U,\lambda _{2} =-\left(\mu +A\right),\left(a_{1} ,a_{2} \right)>0\). Then, all roots of the characteristic equation above have negative real part. Therefore, the endemic equilibrium of Equation (1) is locally asymptotically stable if \(R_{0} >1\).

4. Numerical simulations

In this section, we would show numerical simulations for the proposed Model A. This analysis was carried out in Python with the following initial conditions and parameter values \(S = 600; I = 100; T = 0; R = 0.\) In Table 2 we give the detailed explanation of the parameters, variable and values used in the Model.

Table 2. Table to test captions and labels

Parameters	Description	Value/days	Source
\(\Lambda \)	Recruitment rate	45	Estimated
\(\beta _{1} \)	Effective contact rate	0.00093	[4]
\(\beta _{2}\)	Saturation treatment rate	0.0031	[5]
\(\eta \)	Enhancement factor	0.2	Assumed
\(P\)	Proportion of infected individuals joining either the class R or T	0.04	[5]
\(\gamma \)	Rate of treated individuals from infection	0.7	[15]
\(\sigma \)	Rate at which treated individuals move to recovered class	0.9	[7]
\(\omega \)	Rate at which recovered individuals move to susceptible class	0.8	Assumed
\(\mu \)	The natural death rate in all compartments	0.063	[13]
\(\kappa \)	Educational adjustment	0.012	[4]

4.1. Results

Using Python software, we were able to get the following results. The plots are considered for a deterministic model ‘SITR’, we computed basic reproduction number \(R_{0}\) as an important threshold parameter that has played great role to predict if the diarrhea disease will die out or persist within the population. The detailed explanations of parameters used are under each figure.

Figure 2 shows the numerical simulation solutions of the system model A with different \(R_{0}\) values with the following parameter values \(\Lambda = 45; \beta _{1} = 0:0035; \beta _{2} = 0:05; \eta = 0:02; P = 0:04; \gamma = 0:7; \sigma = 0:9; \omega = 0:8; \mu = 0:063; k = 0:012.\) From figure 2a the infected individuals decrease to zero. This means that, the more saturation treatment the more disease dies out of the population. From Figure 2b as long as we reduce the education and saturation treatment rate, the infected class increases, however when we increase the rate at which individuals are treated to 99%, the infected class goes to zero i.e the disease dies out as time t increases. This shows us that the saturation treatment has an effect to the spread of diarrhea. From figure 2c , we have seen that as long as we increase the education k, then saturation treatment rate; has an effect on treatment as \(\beta _{2} \mathrm{>}\mathrm{>}\beta _{1}\) increasingly. Thus, by varying the treatment, the recovered individuals increase from 19-37 in 1 day. This means that, the more saturation treatment the more people recover, and consequently the disease will die out of the population. From Figure 2d with \(R_{0} \mathrm{>} 1\) means that the diarrhea disease will increase and cannot be eradicated from the population and it is expected that the disease will remain endemic.

5. Conclusion

In this paper, we have analysed the disease-free equilibrium (\(E_{0}\)), the endemic equilibrium (\(E_{1}\)), and the basic reproduction number (\(R_{0}\)). We have shown the local and global stability of the disease-free equilibrium point. We have shown that the endemic equilibrium point is locally asymptotically stable if \(R_{0}>1\) and this means that the diarrhoea disease will persist in population. We have also seen that if \(R_{0}< 1\) the diarrhoea disease will die out of the population after some time. Furthermore, we used python to carry out numerical simulations, and the results show that the maximum possible treatment (saturation treatment rate) affect the dynamics of diarrhea disease. The more saturation treatment, given to individuals, the more the disease dies out of the population. Treatment mention in this project means, giving the infected population the anti-diarrhoea drugs that is: Imodium (loperamide) and Pepto-Bismol or Kaopectate (bismuth subsalicylate). These drugs are helpful in stopping an occasional bout of diarrhoea especially traveler's diarrhoea, which may result from ingesting contaminated food or water while abroad. We recommend the health policy makers that the drugs should be made available to consumers at saturation treatment rate of 99% at a very low cost on time to reduce the dynamical spread of diarrhoea in a community and to avoid epidemics.

Author Contributions

All authors contributed equally to the writing of this paper. All authors read and approved the final manuscript.

Competing Interests

The author(s) do not have any competing interests in the manuscript.