Stochastic and deterministic mathematical model of cholera disease dynamics with direct transmission

In this paper we develop a stochastic mathematical model of cholera disease dynamics by considering direct contact transmission pathway. The model considers four compartments, namely susceptible humans, infectious humans, treated humans, and recovered humans. Firstly, we develop a deterministic mathematical model of cholera. Since the deterministic model does not consider the randomness process or environmental factors, we converted it to a stochastic model. Then, for both types of models, the qualitative behaviors, such as the invariant region, the existence of a positive invariant solution, the two equilibrium points (disease-free and endemic equilibrium), and their stabilities (local as well as global stability) of the model are studied. Moreover, the basic reproduction numbers are obtained for both models and compared. From the comparison, we obtained that the basic reproduction number of the stochastic model is much smaller than that of the deterministic one, which means that the stochastic approach is more realistic. Finally, we performed sensitivity analysis and numerical simulations. The numerical simulation results show that reducing contact rate, improving treatment rate, and environmental sanitation are the most crucial activities to eradicate cholera disease from the community.

Cholera is an acute diarrheal illness caused by the gram-negative bacteria Vibrio cholerae which lives in an aquatic environment (Cui et al. [5]).
The symptom of cholera includes extreme vomiting, dry mouth, irregular heartbeat, painless watery stools, little or no urine output, and low blood pleasure (Fatima et al. [7]). It can be transmitted either by direct and indirect transmission pathways (Fakai [6]). Human-to-human (direct) way of cholera transmission is from the infected individual to the other individuals (touching, biting, and sexual intercourse). Whereas indirect (environment-to-human) way of transmission of cholera is through ingesting vibrio cholera bacteria from contaminated foods and waters (Wang and Modnak [14]). In this paper we concentrate only on the direct transmission path way by incorporating the stochastic nature of the disease.
Mathematical modeling has been an important tool in analyzing the spread and control of infectious diseases and also in making decision as regards the intervention mechanisms for the control of disease (Adewale et al. [1]). Mathematical epidemiology contributed to the understanding of the behavior of infectious diseases, its impacts and possible future predictions about its spreading. Mathematical models are used in comparing, planning, implementing, evaluating, and optimizing various detection, prevention, therapy, and control programs (Bubniakova [3]). Several mathematical models on cholera were developed by different authors. The study done by Codeco [4] focused mainly on endemicity of cholera and suggested two controlling mechanisms. The other study done by Codeco [4] proposed an SIR-B deterministic model by adding an environmental component into the regular SIR model. The study done by Fatima et al. [7] in Nigeria also proposed a deterministic mathematical model for the control of cholera. Additionally, a lot of authors like Fakai [6], Adewale et al. [1], Beryl et al. [2], Javidi and Ahmad [8], and others developed mathematical model of cholera to explore the transmission dynamics and controlling strategies. But none of them consider some stochastic environmental factors that can cause cholera outbreaks like, water, rain fall, air temperature, etc. In this paper we consider environmental factors and develop a stochastic model of cholera dynamics.
The rest of the paper is organized as follows: In Sect. 2, the cholera model is described and formulated in deterministic as well as stochastic approach. In Sect. 3, qualitative analysis and sensitivity analysis of the model are discussed. In Sect. 4, we use MATLAB software to investigate numerical simulation results. Finally, our discussions and conclusions are presented in Sect. 5.

Model formulation and description
The model considers a total human population size (P(t)) and is divided into four compartments, namely susceptible represented by S(t), infected (I(t)), treated (T(t)), and recovered (R(t)) classes. Susceptible individuals (S(t)) are those individuals that cannot be infected, but can get infection sometime in the future, infected individuals (I(t)) are individuals who have developed the symptom of cholera and are able to transmit the disease, treatment class (T(t)) includes those individuals that get treatment at a time t for t > 0 after they have been infected with cholera, and recovered compartment (R(t)) includes those individuals that have recovered from cholera disease and got temporary immunity.
Population in the susceptible compartment will be increased from the recovered compartment with a rate of δ by losing temporary immunity and also with a recruitment rate of π . However, its number decreases by the natural causing death rate of μ and also moving to the infected compartment with the rate of α. Population in the infected compartment will be increased by the contact rate of α, and also its number decreases by the natural causing death rate of μ, cholera causing death rate τ , and moving to the treatment compartment with the treatment rate of σ . Population in the treatment compartment increases from the infected compartment with the treatment rate of σ and decreases with the recovery rate of γ and the natural causing death rate of μ. Population in the recovered compartment also increases by the recovery rate of γ , but its number will decrease by the natural causing death rate of μ and by losing their immunity at the rate of δ.
The model is guided by the following assumptions: the size of homan population is constant, the birth rate and death rate are not equal, all parameters are nonnegative, all individuals are susceptible in the community, therapeutic treatment is applied to the infectious individuals, the treated individuals (individuals that are on treatment) do not transmit cholera disease to the susceptible human population, on recovery there is temporary immunity and there is the disease induced death (disease causing death). By the above descriptions and assumptions, our model is expressed diagrammatically in Fig. 1.
Since the above deterministic model in equation (2) does not consider stochastic environmental factors and lacks realistic conditions, we extended it to a stochastic model.
To extend, we introduce Brownian motion (B i (t)) and the intensity of stochastic environmental factors (β i ) on equation (1) and multiply it by dt. Then we obtain the following  (3): where β 1 , β 2 , β 3 , β 4 ≥ 0 denotes the intensity of Brownian motion and B 1 , B 2 , B 3 , B 4 are independent Brownian motions.

Qualitative analysis
In this section, some basic qualitative behaviors of the model, including the invariant region, positivity of the solution, disease-free equilibrium point, basic reproduction numbers, local and global stabilities of disease-free equilibrium, endemic equilibrium point, stability of endemic equilibrium, and sensitivity analysis, are discussed.

Invariant region
In this subsection we obtain invariant regions of model equations (1) and (2).

Invariant region for deterministic model
To get invariant region, we consider the total population at a time t given by P(t) = S(t) + I(t) + T(t) + R(t). Then differentiating P(t) with respect to time t on both sides, we get Substituting dS(t) dt , dI(t) dt , dT(t) dt , and dR(t) dt from equation (1) into equation (3), we get If there is no infectious individual due to cholera disease, which means (τ = 0, g = 0), then equation (4) becomes By the separation of variables of differential in equality (5), we obtain where C and D are constant. After solving and evaluating equation (6) as t goes to ∞, we have lim t → ∞P(t) = π μ , which implies that 0 ≤ P(t) ≤ π μ . Therefore, the model is biologically meaningful and bounded in the domain = (S, I, T, R) ∈ R 4 + : 0 ≤ P(t) ≤ π μ .

Invariant region for stochastic model Theorem 1
The region is almost surely positive invariant of our stochastic model equation (2).
Proof Let us take K 0 as a large integer so that if (S 0 , I 0 , T 0 , R 0 ) ∈ R 4 + , then every component of (S 0 , I 0 , T 0 , R 0 ) lies in the interval [ 1 k 0 , 1]. For each integer K ≥ k 0 , we can define stopping time: We want to prove P(τ = ∞), which is P(τ < M) = 0 for M > 0, so that it allows us to show lim t→∞ sup P(τ k < 0) = 0. Now consider a Lyapunov function V technically: Then, by using Ito's formula, for Hence, we have Let , then equation (7) becomes Now, rearranging equation (8), we get Then, by integrating both sides from 0 → τ k ∧ M implies that where τ k ∧ M = min{τ n , t}. Taking expectation of the above inequalities yields Since V (X(τ k ∧ M)) > 0, then Now, for τ k , there are some components of Thus, from equation (13) and the above equation, we get From (12)- (14) it follows that Letting k → ∞ and by taking limit sup to equation (16), for all M > 0, we obtain that Therefore, lim t→∞ sup P(τ k < M) = 0, this completes the proof of the theorem.

Positivity of the solutions
In this subsection, we obtain nonnegative solutions for future time with their respective initial conditions. Proof To prove Theorem 2, let us take the first equation of model (1): Then, after solving equation (17) by using the separation of variables and applying the initial conditions, we get After some steps, we get Next, let us take the second equation of model (1): Then, solving equation (19) by using the separation of variables and applying the initial condition, we obtain As t → 0, then e -(μ+d+θ+τ )t → 1, this implies Similarly, we took the third and fourth equations of model (1). Now, taking similar steps to the above gives This completes the proof of the theorem. Therefore, our model solutions of equations (1) and (2) are positive for future time.

Disease-free equilibrium point
To obtain a disease-free equilibrium point, set model equation (1) to zero and there are no infectious individuals in the population, which means I = 0, T = 0, R = 0.

Basic reproduction number
To determine the basic reproduction number (R 0 ) of equations (1) and (2), we use the next generation matrix method by Beryl et al. [2]. We have two basic reproduction numbers (deterministic and stochastic).

Basic reproduction number for deterministic model
In view of that, first let us take the newly infectious class Now, by the principle of next generation matrix approach, we obtain The next step is obtaining the Jacobian matrix of f and v with respect to I at E 0 = ( π μ , 0, 0, 0). Let F and V be the Jacobian matrix of f and v, respectively, then Then, evaluating F and V at a disease-free equilibrium point (E 0 = ( π μ , 0, 0, 0)), we obtain Then FV -1 becomes The eigenvalue of FV -1 can be obtained by Here, by the next generation matrix principle, the largest eigenvalue is the basic reproduction number.
Therefore, our basic reproduction number for the deterministic model is

Basic reproduction number for stochastic model
By taking the infected class of (2), we obtain the basic reproductive number of stochastic approach: Using twice differentiable function of Ito's formula, we can derive our stochastic basic reproduction number. Let us take f (t, I(t)) = ln(I(t)), then its Taylor series expression becomes where ∂f ∂t = 0, ∂f Let a = αS(t)I(t) -(μ + τ + σ )I(t) and b = β 2 I(t), then df t, By applying the chain rule, we get Then df t, By using the next generation matrix, let .
where R D 0 = απ μ(μ+τ +σ ) . Therefore, our basic reproduction number for stochastic model is From equation (28) we see that R S 0 < R D 0 , which means the stochastic approach is more realistic than the deterministic approach.

Local stability of disease-free equilibrium
In this subsection we show the local stability of disease-free equilibrium in the case of deterministic as well as stochastic models. Proof To prove Theorem 3, first we construct a Jacobian matrix of model equation (1).
Then the Jacobian matrix evaluated at E 0 becomes From equation (31) the eigenvalues are evaluated as follows: The characteristic polynomial of equation (32) becomes By factorizing the characteristic polynomial equation, eigenvalues are as follows: Since λ 1 < 0, λ 2 < 0, and λ 3 < 0, but we do not know the sign of λ 4 . However, to be stable all the eigenvalues have to be negative and λ 4 has to be less than zero. This means We know that from equation (24) we have R D 0 = απ μ(μ+τ +σ ) . Hence, equation (34) becomes Therefore, our disease-free equilibrium point is locally asymptotically stable if and only if R D 0 < 1.

Local stability of disease-free equilibrium in the case of stochastic model Theorem 4
If R S 0 < 1, then for any initial values of (S 0 , Proof Let us take F(t, I(t)) = ln I(t), by Ito's formula df t, Integrating equation (35) on both sides, we have Then evaluating equation (36) at E 0 , we obtain where a martingale G(t) = t 0 β 2 dB(t). By the strong law of large numbers of martingale, we have lim t→∞ sup G(t) t = 0 almost surely.
Then divide both sides of equation (37) by t. By letting t → ∞, we get Taking lim t→∞ sup on both sides of equation (38), we obtain Obviously (μ + τ + σ ) > 0, therefore R S 0 -1 less than zero: Therefore, our disease-free equilibrium point is locally asymptotically stable if and only if R S 0 < 1.

Global stability of disease-free equilibrium
Theorem 5 If R D 0 < 1, then E 0 is globally asymptotically stable in .
Therefore, by LaSalle's invariance principle, every solution of model equation (1) with a given initial conditions in approaches to E 0 at t and goes to infinity whenever R D 0 ≤ 0. Hence, disease-free equilibrium is globally asymptotically stable in the region .

Endemic equilibrium point
In this subsection, we obtain the equilibrium point at which no disease is present in the population.
The endemic equilibrium point of our model is denoted by E 1 = (S * , I * , T * R * ). The endemic equilibrium can be obtained by equating all the expressions of model equation (1) to zero: Let us take the second equation of system (44), we get αIS -(μ + τ + σ )I = 0.
Then, solving for S * , we have From the third equation in (44), we solve for I * : Now, solving for R * and T * by substituting equation (45) and (46) into the first equation of (33), we obtain Then, by taking simultaneously the fourth equation of (44) and equation (47), we get Finally, from equation (48), we obtain Then substituting equation (49) into the second equation of (48), we have Now substituting equation (50) into equations (46) and (49), we obtain Therefore, endemic equilibrium point of the model is

Stability of endemic equilibrium point
In this subsection, we discuss the stability of endemic equilibrium E 0 by considering a Lyapunov function L technically as follows: L S * , I * , T * , R * = S -S * -S * ln S + I -I * -I * ln I Now, differentiating both sides of equation (51) with respect to t, we get Then, simplifying equation (52), we obtain the following equation: Now, by simplifying and rearranging (the positive part on one side and the negative part on the other side) equation (53), we have Now, replacing U for the positive terms and V for the negative terms of equation (54), we obtain U = π + αIS * + μS * + αIS + (μ + τ + σ )I * + (μ + γ + g)T * + (μ + δ)R * , and Then, equation (54) is replaced by U and V, and we have If U < V , then dL dt ≤ 0 and dL dt = 0 if and only if our equilibrium points are (S = S * , I = I * , T = T * , R = R * ).
From this, we observe that E 0 is the largest set of compact invariant singletons in {(S = S * , I = I * , T = T * , R = R * ) ∈ : dL dt = 0}. Therefore, E 0 (endemic equilibrium) is globally asymptotically stable in if U < V .

Sensitivity analysis
In this subsection, we obtain sensitivity analysis of the model to determine the effect of each parameter on basic reproduction number (R 0 ). To perform the sensitivity analysis of (R 0 ), we use the following formula: where n i is the parameters of R 0 (Tilahun [13]).

Interpretation of sensitivity analysis
The interpretation of the sensitivity indices given in Table 1 is as follows. Those parameters that have positive sensitivity indices (π, α) have a big contribution to the expansion of cholera disease in the human population if their values are increased by keeping the rest of parameters constant. And those parameters that have negative sensitivity indices (μ, τ , σ )  show a great effect in bringing down the disease from the population if their values are decreased by keeping the rest of parameters constant. Due to the reason that R 0 (basic reproductive number) increases as its parameter value increases, the average number of secondary infection increases in the population; and R 0 decreases as its parameter value decreases, which means that the average number of secondary infection decreases in the human population.

Numerical results and discussion
In this section, some numerical results are presented to demonstrate how change in parameters of the model influences various performance measures of the system. Models (1) and (2) were simulated using MAT LAB software to obtain the following graphs. Additionally, we used S(0) = 12, I(0) = 8, T(0) = 5, and R(0) = 3 as initial values, where the initial values are estimated. In Table 2, the listed parameter values are used for numerical simulation purposes.

Comparison of deterministic and stochastic trends of the model
In Fig. 2, the numerical results of the comparison of deterministic and stochastic trends of the model in the community are displayed by keeping all parameters unchanged. By taking the whole compartments of the model for deterministic figure and by adding a white noise to a deterministic equation, we illustrate a stochastic figure. Now, from this figure, we see that running the simulation results for the stochastic model is slower than for the deterministic model, this is due to environmental factors. The amount of infectious population decreases and the number of individuals who get treatment increases after a certain point of time due to treating the infected people in the community with the rate (σ ) in the stochastic case as well as in the deterministic one. Moreover, the stochastic behavior of the curves shows the real life behavior compared to the deterministic approach. We conclude from this figure that the stochastic solution is closer to the real solution of the cholera model than the deterministic approach. So, using the stochastic model is better than the deterministic one because the stochastic model considers a white noise or stochastic environmental factors.

Effect of contact rate on cholera infected population
In this subsection, we obtain the numerical simulation results of the impact of contact rate (α) on the number of infectious individuals I(t). The simulation results that are displayed in Fig. 3 are obtained by different values of contact rate (α) from α = 0.001 to α = 0.017 and keeping the rest of parameters constant. From Fig. 3, we see that as the contact rate increases the number of infectious individuals in the community increases, which means if the susceptible individuals contact with infected individuals either by shaking hands or eating foods with the same materials this would increase the infectious population.
In addition, we observe that the stochastic curves show a sound wave property due to environmental factors or random behavior of the disease, but in the deterministic model the curve looks like a smooth line, which implies it does not consider any factors in the environment for cholera disease. Hence, the disease persists in the community when there is an increase in the contact rate, even though the rest of the parameters are kept constant. Due to this, healthy workers must control this parameter.

Effect of treatment rate on cholera infected population
We investigated the effects of treatment rate on the number of infectious population. In Fig. 4, the experimental results are obtained by taking different values of treatment rate σ and by keeping the rest of the parameters constant. The simulation results in Fig. 4 reveal We also observe in Fig. 4 that in the case of deterministic approach as well as stochastic approach the decrements of infectious population are obtained as the treatment rate in the population increases. Therefore, increasing the treatment rate σ plays a vital role in the reduction of cholera disease dynamics in the community.

Effects of recovery rate on recovered individuals
In this subsection, we use the impact of recovery rate γ on the size of recovered individuals to display Fig. 5 by varying γ (recovery rate) from γ = 0.03 to γ = 0.7 and keeping the rest of the parameter values fixed. In the case of deterministic approach, Fig. 5 shows that the graph goes up smoothly as the recovery rate γ increases, which means that if the infectious individuals get treatment and recover more with the recovery rate γ , the number of recovered individuals increases in the community. Also, in the stochastic case, it displays that the number of recovered individuals increases as the value of recovery rate increases, with the graph going up and down. These ups and downs show the random behavior of the model. From this, we conclude that the recovered population becomes bigger by increasing the recovery rate (γ ), and the stochastic approach is more advisable than the deterministic one because the stochastic approach considers environmental white noise.

Discussions and conclusions
In Sect. 2, we proposed and described briefly an SITRS cholera disease dynamics model with two approaches: deterministic and stochastic. In Sect. 3, we analyzed the qualitative behavior of the model by obtaining the invariant region, existence of a positive invariant solution set, equilibrium points, the basic reproduction number for deterministic as well as stochastic model; local stability of disease-free equilibrium for both models was checked; the global stability, which is obtained by a Lyapunov function, sensitivity analysis, and their interpretation were studied. In Sect. 3, the two reproduction numbers were obtained by using the next generation matrix method and by using twice differentiable Ito's formula for stochastic reproduction number. Out of these two reproduction numbers the stochastic reproduction number is much smaller than the deterministic one. This implies that the stochastic approach is more realistic (close to the accurate solution) than the deterministic approach, because the stochastic model considers stochastic environmental factors or takes the randomness process. In Sect. 4, the numerical simulation results were discussed and analyzed by using MATLAB software by comparing the deterministic approach with the stochastic one. From our simulation results in Sect. 4, we conclude that increasing cholera treatment rate has a big contribution to eliminating cholera disease from the community, and increasing the recovery rate contributes to the reduction of the infection. Other results that were obtained in this section are as follows: decreasing the contact rate has a big influence on controlling cholera disease dynamics in the community. Therefore, for healthy workers as well as policy makers we recommend improvement of the treatment: by decreasing the contact rate and by increasing the recovery rate we can eradicate the disease from the community.