MARKOV CHAIN MONTE CARLO ANALYSIS OF CHOLERA EPIDEMIC

Several mathematical models have been designed to understand the dynamics of cholera epidemic from which some models considered direct and indirect transmission. In this study a system of ordinary differential equations is developed by splitting the class of infected individuals into symptomatic and asymptomatic infected individuals with the incorporation of water treatment as a control strategy. Theoretically, the developed model is analysed by studying the stability of equilibrium points. The results of the analysis shows that there exist a locally stable disease free equilibrium point, E0 when R0 < 1 and endemic equilibrium, E∗ when R0 > 1. Numerically, the identifiability of parameters is done by least square and Markov chain Monte Carlo methods. Both methods are used as tools to analyze the developed model. The results show that the parameters are identifiable.


Introduction
Cholera is a living testimony of poor sanitary conditions.It is a severe water borne infectious disease caused by the vibrio cholerae bacterium [25].Cholera has short incubation period, from less than one day to five days.It is characterized by severe watery diarrhoea caused by the production of cholera toxin by vibrio cholera bacteria in small intestine and it can cause death within three to four hours if untreated.It is caused by eating food or by drinking unsafe water which are contaminated with vibrio cholerae.It has been proved that pathogenic vibrio cholerae can survive refrigeration and freezing in food supplies [24].
Cholera has been declared by World health organization (WHO) as the public health problem.
Hence, there is a need of finding better ways of dealing with cholera so as to reduce cases of cholera in different countries in the world.According to WHO report, about 1.4 to 4.3 million cases of cholera are reported each year worldwide and more than 140,000 deaths per year are reported due to cholera.The cholera outbreak in Tanzania which began in August 2015 has resulted in over 24,000 cases as of 20 April, 2016 and caused 378 deaths [27].There are several measures which have been suggested by WHO to prevent cholera these measures are environmental sanitation, water treatment, provision of clean water, provision of education on the effect of cholera [26].
A number of studies have been conducted to highlight the spread of infectious diseases in deterministic context.For example, the mechanistic model of the SIRS form, for cholera in [9], explains monthly cholera deaths counts in the twenty-six districts of the former British East Indian province of Bengal during the period 1891-1940.The model incorporated both transmission due to human prevalence via a mass action term and transmission from the environmental reservoir.One of the three models proposed is a two path model which includes a class for severe infectious as well as a class for mild, in apparent infectious.However, the model does not allow for feedback from infected individuals into the environment reservoir.The SIWR model of [19,5] allows for infections from both a water compartment (W) and direct transmission and considers the feedback created by infected individuals contaminating the water.To allow for the possibility of asymptomatic individuals excreting vibrio cholerae to water reservoir.The mathematical model with a compartment for asymptomatic infected individuals developed in [11], considers only direct transmission for disease.
There exist literature that deals with measures of cholera treatment.For example, a mathematical model for the dynamics of cholera with control measures such as educational campaigns, vaccination, sanitation, and treatment as control strategies in limiting disease are explained in [4].The mathematical model considers infected individual as a single compartment.However, it could be better to divide the infected individuals into two groups i.e., asymptomatic and symptomatic infected individuals in order to observe the contribution of vibrio cholerae to the environment from each compartment.The mathematical model (SIR-C) for modeling cholera dynamics with a control strategy in Ghana proposed in [16].The model considers the infected individuals as a single compartment with limited numerical analysis.In this study we extend the deterministic model developed in [16], by splitting infected compartment (I) into two sub groups, symptomatic infected (I s ) and asymptomatic infected (I a ) individuals in order to observe the contribution of vibrio cholerae to the environment from each compartment.

Model Formulation and Theoretical Analysis
We formulate the basic model for the dynamics of cholera with two subpopulation; bacteria (pathogen) and individuals.Individuals are subdivided into four developing compartments S, I s , I a and R, which all of them depend on time t, but the dependency has been dropped for notational convenience.Here S denotes susceptible individuals who contract disease at rate β and the influx of susceptible comes from a constant recruitment rate b and develop to infectious classes at probabilities p, q respectively.The symptomatic infected individuals (I s ) who become new infected from S at a probability p and contribute vibrio cholerae through excretion to the environment at a rate α 1 , dies due to natural death and due to disease at rates µ and d.The asymptomatic infected individuals (I a ) who become new infected from S at a probability q and contribute vibrio cholerae through excretion to the environment at a rate α 2 .R denotes the recovery of I s and I a at rates r 1 and r 2 respectively.The concentration of vibrio cholerae in water is denoted by B. The concentration of bacteria decrease due to mortality rate φ and due to water treatment at a rate δ .The total population is given by N(t) = S(t) + I s (t) + I a (t) + R(t) at any given time t.In formulating the model, the following assumptions are imposed: (1) The population is closed (i.e., there is neither immigration nor emigration), (2) The contribution of both symptomatic infected (I s ) and asymptomatic infected (I a ) individuals to the population of vibrio cholerae in the aquatic environment at the rates α 1 and α 2 respectively, (3) Human birth and death rates occurs at different rates (i.e., b and µ) respectively, (4) Water treatment is included in model as a control strategy, (5) The population is homogeneously mixed i.e., each individual within the population is susceptible to disease.
In this study we assume that, each susceptible individual has equal chance of acquiring cholera through the recruitment rate b and consuming water with vibrio cholerae in the reservoir at the force of infection λ = β B κ + B , where B κ + B is the ratio of vibrio cholerae concetration and κ is the concetration of vibrio cholerae in the water reservoir that will make a possibility of 50% of susceptible population infected.The cholera model can be described by the following deterministic system of nonlinear ordinary differential equations: with initial conditions S(0) > 0, I s (0) ≥ 0, I a (0) ≥ 0, R(0) ≥ 0, B(0) ≥ 0 and p + q = 1.

Computation of SIRB Basic Reproduction Number
In epidemiology a key parameter is the basic reproduction number R 0 , defined as the average number of secondary infectious cases transmitted by a single primary infectious cases introduced into a whole susceptible population [18].To compute R 0 , we use the next generation matrix approach as described in [20].It is obtained by taking the largest (dominant) eigenvalue value (spectral radius) of , where F i is the rate of appearance of new infection in compartment i, V i is the net transition between compartments, E 0 is the disease free equilibrium and X i stand for the terms in which the infection is in progression i.e., I s , I a and B in the model (2.1), that is Using the linearization method, the associated matrix at DFE for F and V can be computed as .
For the inverse of matrix V to exist this condition should hold |V | = 0 and µ, r 1 , r 2 , α 1 , α 2 , δ , µ and d are positive.Hence the next generation matrix becomes , where Then R 0 can now be computed as where

Positivity and Boundedness of Solutions
Model (2.1) can be shown that state variables are non-negative and the solutions remain positive for t ≥ 0. Also, the parameters in our model are assumed to positive.We also, show that the feasible solutions are bounded in a region such that Φ = (S, I s , I a , R, B) ∈ R 5 + .Lemma 2.1.Let the initial values of the parameters be {S(0) ≥ 0, I s (0) ≥ 0, I a (0) ≥ 0, R(0) ≥ 0, where b = bN.
We have that, by separating the variables of Equation (2.3) and integrating, we obtain Hence, By considering the second equation in Equation (2.1) We have that By separating the variables of Equation (2.4) and integrating, we obtain Hence, By considering the third equation in Equation (2.1) We have that By separating the variables of Equation (2.5) and integrating, we obtain Hence, The same method can be applied to the remaining equations in fourth and fifth in Equation (2.1) to obtain Hence, Hence, Therefore, the solution of model system (2.1) is always positive.This completes the proof.
Lemma 2.2.The solutions for the model system (2.1) are contained and remain in the region Φ for all time t ≥ 0 Proof.Consider the total population The integration factor I.F= e µt .The solution becomes where However, for the bacteria variable the boundedness is shown as follows By integrating this equation (I.F= e (δ +φ )t ), we get this solution where A is a constant.Then which implies that N and all other variable (S, I s , I a , R and B) is bounded and all the solutions starting in Φ approach, enter or stay in Φ.This completes the proof.

Local Stability of the Disease Free Equilibrium
Local stability of the DFE can be analyzed using R 0 as the bifurcation parameter, that is locally asymptotically stable if R 0 < 1 and unstable when R 0 > 1.The DFE of the model system (2.1) is given by The DFE E 0 of the system (2.1) is locally asymptotically stable if R 0 < 1 and Proof.The Jacobian matrix of the system at an arbitrary equilibrium is defined by .
The characteristic equation of the Jacobian matrix is given as follows Since the first and fourth columns of matrix above contains diagonal terms, they form eigenvalues λ 1 = −µ which is negative.The other eigenvalues are obtained after reducing the first and fourth columns and their corresponding rows, leading to where The simplification of Equation (2.7), leads to ) and Equation (2.8), can be written as where We can now write Equation (2.9), in the form where To ensure that the remaining eigenvalues of Equation (2.10) have negative real parts, we employ Routh-Hurwiz stability criterion [17].The conditions are M 1 > 0, M 3 > 0 and also, M 1 M 2 > M 3 .So, M 1 is already non-negative and M 3 is non-negative if and only if R 0 < 1.
The disease free equilibrium is locally asymptotically stable if R 0 < 1.We re-write M 3 as when R 0 < 1, M 3 is positive and R 0 > 1, M 3 is negative, under the condition that L is always positive, which is true since the values of the parameters are all positive.This proofs Theorem 2.3 i.e. the disease free equilibrium of the system is locally asymptotically stable if R 0 < 1 and unstable if R 0 > 1.This completes the proof.
Theorem 2.2.The DFE is globally asymptotically stable if R 0 < 1 Proof.For globally asymptotically stable, we use a concept of Metzler matrices proposed by [2].The system must be written as where X ∈ R m denotes the uninfected compartments and Z ∈ R n denotes the infected compartments.The DFE is globally asymptotically stable for the system provided that R 0 < 1 and conditions stated below should holds, H 1 : The system is defined on a positively invariant set Ω of the non-negative orthant.
That is dX dt = M(X, 0), X * is globally asymptotically stable.
Here G = ∂ N ∂ Z (X * , 0) is an M-matrix (the off diagonal element of G are non-negative).
Therefore X = (S, R) and Z = (I s , I a , B) At the disease free equilibrium, Z = B, I s , I a , but Z = 0, so , B = 0, I a = 0, I s = 0, and hence leads to Ṡ = bN − µS and Ṙ = −µR, and their corresponding solutions becomes It is true that R(t) → 0 and S(t) → bN µ as t → ∞, regardless of the values of R(0) and S(0).
Hence, we conclude that the system is globally asymptotically stable at the equilibrium point bN µ , 0 , thus H 1 are satisfied.Again the matrix N(X, Z) is given by .
Hence, we observe that N(X, Z) is less than zero i.e., N(X, Z) < 0. Therefore, H 2 is not satisfied.
We conclude that, the disease free equilibrium may not be globally asymptotically stable.This completes the proof.

Endemic Equilibrium Point and Local Stability
The endemic equilibrium points of the model system (2.1) is given by E * = (S * , I * s , I * a , R * , B * ) with I s = 0, I a = 0 and B = 0.It can be obtained by equating the RHS of each equation of the model system (2.1) equal to zero, which exists for R 0 > 1.
Theorem 2.3.If R 0 > 1, the endemic equilibrium E * exists and is locally asymptotically stable Proof.The local stability of the endemic equilibrium is established from the eigenvalues of the Jacobian matrix evaluated at endemic equilibrium, that is From which it is observed that λ = −µ < 0 is an eigenvalues.The other four eigenvalues can be obtained from the characteristic polynomial of the 4 × 4 block matrix.
The remaining condition is to have a stable system, that is Det(J 1 ) > 0.
The easier way to find the determinant of J 1 is to expand using column 2. Hence, we have where For Det(J 1 ) > 0, this condition should hold D 1 > D 3 .
Since tr(J 1 ) < 0 and Det(J 1 ) > 0. We conclude that the system is stable.Thus implying that, from Equation (2.2), R 0 > 1.Therefore, the endemic equilibrium is locally asymptotically stable when R 0 > 1.This completes the proof.where

Global Stability of Endemic Equilibrium Point
By direct substitution to Equation (2.11) we get (2.12) Also, the endemic equilibrium of the model system (2.1) is given by E * = (S * , I * s , I * a , R * , B * ).It can be obtained by equating the right hand side of each equation of the model system (2.1) equal to zero.Thus (2.13) Therefore, Equation (2.12) deduce to Putting the positive and negative terms together in the system (2.15) we obtained where Then, from Equation (2.16), if M < Z, then dV dt will be negative, implying that dV dt < 0. However, it follows that Therefore, the largest compact invariant set in is the singleton E * , where E * is the endemic equilibrium of the model system (3.1).By LaSalles's invariant principle, it implies that endemic equilibrium E * is globally asymptotically stable in Ω if M < Z.This completes the proof.

Numerical Results and Discussions
In this section, simulation and parameter estimation of the developed model in (2.1) is carried using least square and adaptive Markov chain Monte Carlo methods as in [13].Data are created by solving ODEs and then corrupted it with relative Gaussian noise whose standard deviation is 0.5.The parameter values used are literature values and by substituting these values to (2.1), the simulated ODEs (2.1) leads to the results shown in Figure 2. From Figure 1, we see that, susceptible variable is decreasing this is due to the fact that some of its members are immigrating to I a and I s compartments.As the time goes both I a and I s are increasing and later, decrease after a period of time this is due to control measures taken, R is increasing exponentially this implies that all individuals reaching the compartment R will never come back to the system and are supposed to remain within it.The parameters described in model system (2.1) were estimated by least-squares method.
This involves minimizing the sum of squares of residuals.The results are shown in Table 1.
From Table 1, it is observed that estimates of the parameters are indeed close to the literature Estimates 0.000066 0.32 887729 1.4 0.092 0.000586 0.045 0.022 0.12 0.46 0.644 0.27 0.0000402 values which, in this case, are treated as true values.This implies that, least square method performs well to the cholera model developed.
In MCMC parameter sampling, to know weather our chains have converged or not, we use assessment methods such as summary of MCMC, trace plots, scater plots, marginal posterior distribution and autocorrelation functions.The initial values used are S = 1000000, I s = 1, I a = 1, R = 1, B = 1 and we generated 100000 samples using initial covariance of 0.1.Also, we calculated the basic reproduction number R 0 using true values, estimated values and MCMC mean.The R 0 value are 5.7, 6.5 and 5.01 respectively.

A Summary of an MCMC Object
The results gives the posterior means, standard deviations and posterior quantiles for each chain and convergence diagnostic.From Table 2, the mean values are close to the least square  From Figures 2 to 4, we observe that the chains seem to be stationary, there is no obvious trend or stuck, which implies good mixing of chain.

Autocorrelation Function
The autocorrelation functions measure how well the MCMC sampler performs by measuring the autocorrelation between parameters θ i and θ i+q at lag q.The smaller the autocorrelation values, the better mixing of the chains.The autocorrelation values for Figures 5 and 6 are decreasing exponentially and stabilizing around zero, this proves that the parameters are identifiable.From the MCMC figures, we get the information related to correlation, uncertainty, identifiability of parameters, convergence of Markov chain to the target distribution etc [23].The distributions that are skewed to the left have a negative coefficient of skewness and that skewed to the right have a positive value and the skewness for normal distribution is 0 and the kurtosis for normal distribution is 3 [7].A ratio greater than 3 indicates more values in the neighborhood of the mean and a ratio less than 3, indicates that the curve is flatter than the normal.See, Table 4 and Figure 7 for more details.However, it should be noted that parameters for ODEs can take asymptotic properties of any distribution.

Predictive MCMC Plots
We check the accuracy of the model through prediction plots.From Figure 8, the model predicted the data at 95% posterior limits which seen with the gray colour around the model solution.The variance of predictive distribution reflects the predictive accuracy of the model.

Conclusion
In this paper, we formulated a new SIRB epidemic model by splitting the infected compartment into two classes (I s and I a ) with the aim of modeling cholera epidemics.We have derived The basic reproduction number (R 0 ) is calculated and its value is seen to be 5.7 which agrees with the one found in [22] and is greater than 1.We observe that the disease was capable to invade susceptible population, but its spread was possibly stopped due to intervention and other measures taken.
The model parameters have been estimated by applying least squares estimation with the aim of fitting the SIRB ordinary differential equations to data.Also, Markov chain Monte Carlo (MCMC) method is used to estimate unknown parameters and other characteristics of the target posteriors by generating samples.Graphical representations are presented to illustrate and support the analytical results.The predictive distributions generated predicted the model to a large degree of accuracy.Finally, it is observed that both least squares and MCMC methods performed well to the cholera model developed.

Theorem 2 . 4 .
If R 0 > 1, the endemic equilibrium E * of the model system (2.1) is globally asymptotically stable.Proof.To establish the global stability of endemic equilibrium E * , we construct the derivative of positive Lyapunov function V as follows;

FIGURE 1 .
FIGURE 1.Time evolution of the susceptible, symptomatic, asymptomatic infected, recovered and bacteria classes.

FIGURE 2 .
FIGURE 2. Trace plots of estimated unknown parameters b and β using MCMC method.

FIGURE 3 .FIGURE 4 .
FIGURE 3. Trace plots of estimated unknown parameters p and q using MCMC method.

FIGURE 5 .
FIGURE 5. plots of Autocorrelation function of b and φ with 150 lags.

FIGURE 6 .
FIGURE 6. plots of Autocorrelation function of δ and d with 150 lags.

TABLE 4 . 1 FIGURE 7 .
FIGURE 7. Posterior distribution of unknown parameters with 10000 iterations by defining the normal with large variance as the prior distribution on the shape parameter.
Substituting the derivatives from Equation (2.1) to Equation (2.6) we get

TABLE 1 .
Estimated cholera epidemic model parameters by least square method.

TABLE 2 .
A summary of an MCMC of parameter value.
[6]h iteration[6].Through this we check whether the chains get stuck in a certain areas of the parameter space, which indicates bad mixing.