Dynamics of a Stochastic and Deterministic SVIQRS Cholera Epidemic Model

In the present research paper, deterministic and the corresponding stochastic mathematical models describing the dynamics of cholera epidemic are presented by incorporating vaccination. The total population size of the model is divided into ﬁve compartments namely Suscep-tible, Vaccinated, Infected, Quarantined for treatment and Recovered class. Initially, the cholera model is developed, and is determined by a deterministic approach. Since this deterministic approach is not considering either environmental factors or the randomness process of the dynamics, a corresponding stochastic approach has been introduced. The model equations of both deterministic and stochastic cases have been proved to be positive and also bounded. Furthermore, for both the models, mathematical formulations of the basic reproduction numbers are developed by employing the next generation matrix method. The analysis shows that the basic reproduction number for deterministic approach is much greater than that for the stochastic one. Finally, numerical simulations are also performed. The simulation study has revealed that combination of a decrease in contact between infected and susceptible individuals, increasing vaccination coverage, creating awareness to reduce contact rate, increasing recovery rate with proper treatment, and environmental sanitation are the most basic control strategies so as to eliminate cholera disease from the community.


Introduction
Cholera is a severe water and food borne infectious disease.Cholera is caused by Vibrio cholera bacterium, which lives in an aquatic environment.The bacterium transmits cholera directly from human to human and also indirectly from environment to human [1].It is true that an individual infected with the disease may show or may not show symptoms.Some of its symptoms are watery diarrhea, vomiting and leg cramps.If an infected individual is not treated timely then he will suffer with acidosis, dehydrated and circulatory collapse.This situation may lead to a death within a time period of 12 to 24 hours [2] .Some studies and experiments proved that a recovered individual is immune to the disease for a period of 3 to 10 years.Recent researches recommend that immunity can be lost after a period of weeks to months [3].Thus, dynamics of cholera is caused due to multiple interactions among human hosts, pathogens, and environment [4].Between 2007 and 2019, several cholera outbreaks occurred in various countries including Angola, Haiti, Somalia, Congo, Zimbabwe and Yemen [5].Recently, the cholera has broken out in Yemen in a big way [6].The epidemic started in October 2016 and ended by Mach 2017.However, on 27th April 2017 the epidemic returned.Problems in infrastructure, water and sanitation facilities allowed the disease to spread Invalid source specified.In 2018, 2990 deaths and 499 447 positive cases of cholera were registered.According to [7], though the outbreaks are going on in various countries there has been a significant down trend in case of cholera.During 25 April and 6 June 2019, as many as 424 cholera cases and at least 15 deaths have been reported in Ethiopia.The most affected region in Ethiopia is Amhara with 198 cases, followed by Oromia with168 cases, Somali with 33 cases, Addis Ababa with 15 cases and Tigray with 10 cases.Of these, 13 cases were caused due to cultural activities: 5 in Oromia, 4 in Addis Ababa, 2 in Amhara and 2 in Tigray [7].
The disease cholera has a huge impact on public health, economy and also social development.Cholera has been an important subject of research studies in experimental, clinical and theoretical areas.Generally, it is believed that some strategies like quality water, sanitation, and hygiene will control cholera.But, there have been numerous examples of its existence and transmission in spite of application of the afore mentioned control strategies.Therefore, vaccination has been recommended as an additional strategy to control cholera.However, oral cholera vaccine (OCV) is more effective than early-generation parenteral vaccines to control cholera [8].The natural immunity that resists cholera infection to effect individuals depends on vaccination and some other factors [9] So far, a good number of mathematical models have been developed to understand the dynamics of cholera.The models are mathematically analyzed too over the time [ [10], [11], [12], [13], [14]].
In [10], a SIR model is introduced dividing whole human population in to three compartments viz., Susceptible, Infectious, and Recovered and dividing bacterial concentrations in to two compartments viz., hyper infectious and less-infectious.However, the infectious individuals are further subdivided into two compartments as the individuals are asymptomatic and symptomatic.In [11] yet another SIR type of model is presented including some control strategies viz., public health educational campaigns, vaccination, quarantine and treatment.These strategies were shown to be very promising in controlling the spread of the disease.In [12], SIR-B type model is introduced in which vaccination and disinfection as control measures of cholera have been included.In this research it is pointed out that the immigration is a potential factor for the propagation of the disease.In all these models afore mentioned, it is well proved that the cholera outbreak can be kept under control if both vaccination and treatment are administered effectively in the population.In [13] is was shown that by strengthening both prevention and control strategies sufficiently the epidemic can be eliminated from the populations.In [14] the authors studied stochastic and deterministic mathematical models of cholera disease and described the dynamics with direct transmission.They showed that the disease-free equilibrium point E 0 is globally asymptotically stable whenever R 0 < 1, and as a result the disease dies out.Finally, they showed that a decrease in the contact rate has huge influence on cholera disease so as to eliminate it from the community.Considering all the afore mentioned models and their features here a new five compartmental model SV IQRS is proposed.This model is an extension of that given in [14].Here the considered extensions are (i) vaccination and (ii) recovered individuals become susceptible again due to loss of immunity.Deterministic and the corresponding stochastic versions were considered and were compared.The details of the models such as assumptions, mathematical analysis, simulation studies etc. have been incorporated in the subsequent sections.

Description and formulation of modified model 2.1 Model Assumption
The total human population of the model at any time is divided into five groups based on their disease status: Susceptible S(t) , Vaccinated V (t) , Infected I(t) , Quarantined for treatment Q(t) and Recovered R(t).The susceptible class S(t), consists of individuals of all age groups of the population who have not come into effective contact with the Vibrio cholera.The vaccinated class V (t), consists of individuals who had been vaccinated and still possess partially immunity against cholera.The infected class I(t), consists of individuals who are infected with cholera and are capable of propagating the disease to susceptible humans.The quarantine class Q(t), consists of infected individuals.These infected individuals are subjected to stay in quarantine for a specified period of time.In quarantine the infected individuals are isolated and are provided medication.The recovered class R(t), consists of individuals who are successfully treated.These individuals are sufficiently immune against the disease.Population of S(t) class increases with a constant birth rate or recruitment rate Π.Humans progress into S(t) from R(t) with a rate η as they lose the immunity that was acquired from the treatment.However, the population of S(t) decreases due to vaccination and also infection.The class V (t) is increased with the rate of vaccination ϕ i.e., with this rate humans progress into V (t) from S(t).However, some population of V (t) are imperfectly vaccinated and they will go back to S(t) a rate θ.Here 0 ≤ κ ≤ 1 is the protection efficiency of the vaccination and 1 − κ is the risk of infection due to vaccination inefficiency.Thus, the vaccination protection efficacy is considered as 100 percent if κ = 0 and 0 percent if κ = 1.Due to the fact that vaccines are not fully immune for cholera epidemic disease, therefore, the vaccines can still become infected again.As a matter of fact, the vaccination against cholera disease does not stop completely being infected but it reduces the probability of being infected.Hence, we assumed the law of mass action interaction between vaccinated and infected compartments [15] [16].Thus, humans enter into I(t) from S(t) at a rate βτ and from V (t) at a rate βτ κ.Here, τ is the contact rate of susceptible individual with infected, β is the probability that a contact results in the infection and thus βτ is the infection rate or rate of propagation of infection.However, individuals of I(t) progress to Quarantine Q(t) at a rate δ for isolation, treatment and medication.Humans are recruited in to Quarantine compartment Q(t) from infected class I(t) at a rate δ.However from Q(t) , after getting successful treatment and losing the infection, individuals will move to the recovered class R(t) with a rate φ.Also, the recovered individuals will lose the immunity after some time and they become susceptible.In all the classes the populations are expected to decrease due to natural death with a rate µ.Additionally, population sizes of I(t) and Q(t) will decrease at the disease-induced death rates of ξ and ψ respectively.Considering the definitions, assumptions, and inter-relationships among variables and parameters, the basic dynamics of cholera is illustrated in the form of a flow diagram as shown in 1.
Based on the model assumption and the Schematic diagram the model equation is formulated with initial condition: The deterministic model of equations 1 does not consider the effect of randomly fluctuating environment or environmental factors, including contamination of food and water, improper sanitation, malnutrition, occupational risk.To include these environmental factors, it is appropriate to incorporate white noise in each of the model equations 1.Now, let some stochastic environmental factor acts simultaneously on each individual in the population.Let W i (t) be the mutually independent standard Brownian motion with W i (0) = 0 and let i where i = (1, 2, 3, 4, 5), are the intensities of white noise.By introducing these stochastic perturbations, the stochastic version corresponding to the deterministic model given in 1 takes the following form: [17,18]: (2) 3 Qualitative Analysis of the Model equations

Feasible region
Here in this section, the invariant region in which the solutions of the system of equations given in 1 are bounded will be obtained.Now, on differentiating the total population N (t) = S(t) + V (t) + I(t) + Q(t) + R(t) with respect to time t and substituting into Eq. 1, the simplified equation can be obtained as : Initially there is either no infection or that is negligible I ≥ 0 and also the disease induced death rate satisfies ξ ≥ 0. Thus, without loss of generality Eq. 3 can be re-expressed as dN (t) dt ≤ Π − µN which on solving using variables separable method gives: Further, it can be observed that N (t) → Π µ as t → ∞.That is, the total population N (t) takes off from the initial value N (0) at the beginning and ends up with the bounded value ( Π µ ) as time grows to finitely large.Thus, it can be concluded that N (t) is bounded i.e., 0 ≤ N (t) ≤ ( Π µ ).Thus, the solution set of the system of model equations 1 enters and remains in the feasible region: Therefore, the system of model equations 1 is biologically well posed and mathematically meaningful.Hence, it is appropriate and sufficient to study the dynamics of the model variables in the invariant region Ω.

Non-negativity of the solution:
It is assumed that the initial conditions or values of the model variables are nonnegative.Also, it has been shown that solutions of the model equations are positive.
Then, solutions of the system of model equations 1 are positive.
Proof: Positivity of S(t): Consider the first differential equation of 1 as dS dt = Π + θV + ηR − µS − βτ SI and that without loss of generality can be expressed as dS dt ≥ −(µ + βτ I)S.Now, following the method of separation of variables it solved to obtain its solution as S(t) ≥ Ae −(βτ +µ)t where A = S(0) > 0 is an integral constant.It can be observed that S(t) > 0 as t goes to ∞.Thus, it can be concluded that S(t) is a positive quantity.Positivity of V(t): Consider the second equation of 1 as dV dt = ϕS − (µ + θ)V − βτ κV I and that without loss of generality can be expressed as dV dt ≥ −(µ + θ + βτ κI)V .Now, following the method of separation of variables it solved to obtain its solution as S(t) ≥ Ae −(βτ κ+θ+µ)t where A = V (0) > 0 is an integral constant.It can be seen that V (t) > 0 as t goes to ∞.Thus, it can be concluded that V (t) is a positive quantity.Similarly, by applying the same technique for remaining equations we obtained: Thus, the system of model equations 1 are positive for all t ≥ 0. Therefore, the model is meaning full and well posed in Ω.

Disease Free equilibrium (DFE)
In order to find the disease free equilibrium DFE point of the model, the right hand sides of the system of equations 1 are equated to zero.Then the resultant equations are evaluated at I = Q = R = 0 and they are solved for non-infected and non-carrier state variables.Thus, coordinates of disease free equilibrium point are obtained as: , 0, 0, 0 .
3.4 Threshold quantity (R 0 ) 1. Threshold quantity for deterministic model Here, the threshold parameter that governs the spread of disease known as the basic reproduction number is obtained.It is nothing but the spectral radius of the next-generation matrix [19].For the purpose the system of model equations 1 are rearranged starting with those representing newly infective classes.
Now, using the principle of next-generation matrix method f i and v i are obtained as: Now partially differentiating the variables f i and v i with respect to I and evaluating at the disease free equilibrium point and then the substitution of S = (θ+µ)Π µ(θ+ϕ+µ) , and V = ϕΠ µ(θ+ϕ+µ) reduces the Jacobian matrices to: Now, the product of the matrices F and V −1 can be computed as: Now it is easy to identify eigenvalues so as to determine the required basic reproduction number R 0 .It is nothing but the spectral radius or largest eigenvalue of the matrix F V −1 .Thus, the eigenvalues are computed by evaluating det(F V −1 − λI) or equivalently solving the following matrix equation: However, it is straight forward to identify the largest eigenvalue here as λ = (θ+µ)βτ Π+βκτ ϕΠ µ(θ+ϕ+µ)(δ+µ+ξ) .It is the spectral radius or the threshold value or the basic reproductive number.Thus, it can be conclude that the basic (effective) reproduction number of the deterministic model is given by: 2. Threshold quantity for stochastic model By taking the infected class of the system of model equations 2 as: Using Ito's formula for twice differentiable function f (t, I(t)) = ln(t, I(t)), its expansion in Taylor series is: Where Then equation 8 becomes: Let h 1 = βτ I(t)S(t) + βκτ I(t)V (t) − (δ + µ + ξ)I(t) and h 2 = β 3 I.
Now, by using chain rule, it is obtained as: dtdt = 0, dtdB(t) = 0 and Thus Jacobean matrix is obtained as: Then, disease free equilibrium of the model system 2 at S = (θ+µ)Π µ(θ+ϕ+µ) and V = ϕΠ µ(θ+ϕ+µ) are simplified to: The product of the matrices F and V −1 can be computed as: The Eigen value of By the principle of next generation matrix, the dominant Eigen value is named as the basic reproduction number.Hence, the basic reproduction number for the stochastic case is given by: Similarly, the relation between R S 0 and R D 0 is found as: It is because the stochastic version is much closer to reality than the deterministic one.

Local Stability of Disease Free Equilibrium
Theorem 2.: Disease free equilibrium E 0 of system of deterministic model equations 1 is locally asymptotically stable if R 0 < 1 and unstable if R 0 > 1. Proof: Now, the Jacobian matrix of the deterministic model equations 1 at the disease free equilibrium E 0 reduces to the form as : Here and Now, to find eigenvalues of J(E 0 ) let the characteristics equation det(J(E 0 ) − ψI) = 0 be expanded and be simplified as : (11) From the Jacobian matrix of 11 we obtained a characteristic polynomial: Where a = (θ+µ)βτ Π+βκτ ϕΠ µ(ϕ+k2) As required by the principle of Routh-Hurwitz criteria, equation 12 will have negative real roots if and only if the following conditions hold true: 3 > 0 .Thus, DFE E 0 of the deterministic system of the differential equations 1 is locally asymptotically stable if R D 0 < 1 and unstable if R D 0 > 1. Theorem 3. Disease free equilibrium E 0 of system of stochastic model equations 2 is locally asymptotically stable if lim sup Now, substitution on k 3 = (δ + µ + ξ) reduces df (t, I(t)) to: Now, on integrating both sides of equation 13 gives: Equivalently, it can be expressed in terms of inequality as: Here in 14 the martingale is given by By applying the strong law of martingale, almost surely it is obtained as lim sup t→∞ G(t) t = 0 .Then, dividing both sides of equation 14 by t and evaluating the limit as t → ∞, it is obtained as: Therefore, disease-free equilibrium point is locally asymptotically stable if and only if R S 0 < 1.

Endemic Equilibrium
The endemic equilibrium points E * = S * , V * , I * , Q * , R * are defined as the steady state solutions whenever the disease persists in the population.Mathematically, the endemic equilibrium points can be obtained by setting rates of changes of variables with respect to time, given in model equations 1, to zero.That is, first set dS dt = dV dt = dI dt = dQ dt = dR dt = 0 , then solve the resultant model equations for state variables in terms of the parameters, and after performing some algebraic operations finally obtain the equilibrium as follows:

Sensitivity Analysis
Here, the sensitivity analysis is performed by using the method given by [20].This analysis helps to identify those parameters which show high impact on the disease propagation in the community.Sensitivity indices of the parameters are obtained by using the normalized formula Υ R0 mi = ∂R0 ∂mi * mi R0 , where m i is any parameter in R 0 and are tabulated in Table 1.
For m = β, For m = θ, The physical interpretations of sensitivity indices given in 1 can be made as  the average number of secondary infections increase in the population.In the similar lines, R 0 decreases as its parameter value decreases, which means that the average number of secondary infection decreases in the human population.

Discussion
In this section, we explore the comparison between deterministic and stochastic approach for cholera model in terms of the effect of probability of contact rate, treatment rate and vaccination rate on the infected human population.
Similarly, the effect of treatment on the quarantine population and the effect of recovery on the recovered human population are presented.In Figure 3 , the comparative results of both deterministic and stochastic trends of the disease in the community are displayed, while all parameters are kept constant or unchanged.Population sizes of all the compartments of the deterministic model are shown in Figure 3(a).By adding white noise to the deterministic model equations, the corresponding stochastic model is obtained and its population sizes are illustrated in Figure 3(b).From Figures 3(a) and 3(b), it can be seen that the simulated populations for the stochastic model run slower than that for the deterministic model, which is because of the environmental factors.Further it can be observed from these simulated figures that after certain point of time the infectious population decreases while the treatment population increases, since the treatment is introduced with a rate δ both in the stochastic and deterministic cases.Moreover, the simulated curves show that the stochastic approach is much closer to the real life behavior when compared with the deterministic one.It can be concluded from these figures that the stochastic solution is closer to the real solution for the cholera model than the deterministic approach.So, it can be suggested that the usage of stochastic model is better than that of the deterministic one because the former considers a white noise or stochastic environmental factors.In this section, the impact that the probability of the contact rate parameter β on the number of infected individuals variable I(t) is investigated.In Figures 4(a) and 4(b), the numerical results obtained by varying the value of β while keeping other parameters fixed are presented.In case of deterministic model given in Figure 4(a), whenever the value of β is increased from 0.1 to 0.2, there is a significant and regular increase in the number of infected individuals.Moreover, when β attains the value 0.3, the size of infected population quickly increases and manages to stay higher than in the foregoing two cases.On the other hand, as depicted in Figure 4(b), these results of the stochastic model are also increasing and are also maintaining their zigzagging property due to the randomness behavior.However, the overall outcome is that the number of infected people is still increasing significantly with increase of the value β .Therefore, it can be concluded that even if other parameters are kept constant, the disease expands in the community whenever there is an increase in the probability of contact rate.In Figure 5, the numerical results obtained, by varying the value of contact rate τ , are plotted.The parameter τ is varied and the other parameters are fixed so as to investigate impact of τ on the number of infected population I(t).It is easy to observe that, in both deterministic 5(a) and stochastic 5(b) cases, the curve runs higher for larger values of τ implying that the disease expands in the community.That is, the number of infected people increases together with increasing rate of contact.Here, the effect of treatment rate on the number of infectious population is investigated.In Figure 6, the experimental results are obtained by assigning different values for treatment rate δ while keeping the remaining parameters at constant.The simulation results reveal that an increase in the value of treatment rate leads to decrease in the population of infectious individuals.
In other words, in case of both deterministic and stochastic approaches, the infectious populations are decreased as the treatment rates in the population are increased.Therefore, increasing the treatment rate δ plays a vital role in the reduction of cholera disease dynamics from the community.Here, impact of the recovery rate φ on the population size of recovered individuals is analyzed.In Figure 7, the sizes of recovered individuals are display by varying the recovery rate φ from 0.1 to 0.4 and while keeping the other parameters fixed.In the case of deterministic approach displayed in Figure 7(a), it can be observed 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 will increase 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 reflect the random behavior of the model.From this, it can be concluded that the recovered population becomes bigger by increasing the recovery rate φ.Also, the stochastic approach is more advisable than the deterministic one because the former approach considers environmental white noise.In Figure 8, numerical results obtained by varying the values of vaccination rate parameter ϕ and fixing the other parameters are presented.In the deterministic case given in Figure 8(a), the numbers of infected individuals decrease as the time progresses.The curves of I(t) are slightly higher for smaller values of the vaccination rate ϕ which is obviously expected.Therefore, vaccination of the target population has a significant contribution in eliminating the disease.On the contrary, the stochastic case exhibits a very different behavior and the same is displayed in Figure 8(b).For the first few moments, the number of infected people I(t) is seen to grow upwards for all values of vaccination rate.However, the value of I(t) turns to go down as time goes by.Here, it is worth noticing that the increment in the value of the vaccination rate plays a vital role in eliminating the disease in the stochastic case as well.

Conclusion
In this paper, an SV IQR mathematical model has been developed so as to describe the dynamics of cholera diseases in two different approaches viz., stochastic and deterministic.In section 3, Qualitative behaviors of both these models are analyzed which include the invariant region and positivity of the solution set.Further, it is shown that the equilibrium points exist.Also, the basic reproduction number is constructed in terms of parameters.Local stability of disease-free equilibrium for both these models is verified.Also, the global stability is verified using the Lyapunov function technique.Sensitivity analysis of the models is conducted and useful interpretations were drawn.Two reproduction numbers were obtained by using the next generation matrix method and by using twice differentiable Itos formula for stochastic reproduction number.Of these two reproduction numbers the stochastic one is much smaller than the deterministic one.This implies that the stochastic approach is more realistic or close to the accurate solution than the deterministic approach, because the former model considers stochastic environmental factors or takes the randomness process.In section 4, the sensitivity analysis on the basic parameters is conducted and the physical interpretations of the analysis are presented.In Section 5, numerical simulations are carried out using MAT LAB software and the results of stochastic and deterministic approaches are compared.The paper ends in Section 6, with the inclusion of useful conclusions drawn from the numerical simulations.These results show that the number of infected people keeps decreasing if one carefully combines vaccination with appropriate treatment.It is also observed that the stochastic model mimics and represents better the real-life phenomena when compared to deterministic approach.Moreover, the impact of the parameters β, τ, φ and ϕ is also investigated in case of both the deterministic and stochastic models.
It is observed that increasing the contact rate τ and probability of contact rate β contributes to the spread of the disease whereas decreasing contact rate τ , increasing vaccination rate ϕ and increasing the value of the recovery rate φ will play vital roles in eliminating the disease from the community.Finally, it is clear that the transmission of cholera depends on the probability of contacts and hence consideration of random effects in developing these models makes the transmission dynamics of cholera disease more realistic.

Fig. 1
Fig. 1 Compartmental structure and flow diagram of the model

Theorem 4 .
If R 0 D < 1, then the DFE is globally asymptotically stable in the feasible region Ω.Proof: Now let the Lyapunov function be constructed technically as L(t) = Π + βδ + ϕ µ I(t).The differentiation of L(t) with respect to time and substitution of dI dt = [βτ S + βτ κV − k 3 ] I give rise to dL dt = Π + βδ + ϕ µ [βτ S + βτ κV − k 3 ] I.This equation with our loss of generality, by replacing S and V with the corresponding initial values, can be expressed as an inequality as dL dt
Deterministic case (b) Stochastic case

Fig. 4
Fig. 4 Effect of probability of contact rate β on I(t) for deterministic 4(a) and stochastic 4(b) model.

Fig. 5
Fig. 5 Impact of contact rate τ on I(t) for deterministic 5(a) and stochastic 5(b) model.
Deterministic case (b) Stochastic case

Fig. 6
Fig. 6 Effect of recovery rate φ on I(t) for deterministic 6(a) and stochastic 6(b) model.
Deterministic case (b) Stochastic case

Fig. 7
Fig. 7 Effect of recovery rate φ on I(t) for deterministic 7(a) and stochastic 7(b) model.

Fig. 8
Fig. 8 Effect of vaccination rate ϕ on I(t) for deterministic 8(a) and stochastic 8(b) model.

Table 1
Sensitivity index table

Table 2
Parameter Values for Cholera Model