Continuous Time Markov Chain Model for Cholera Epidemic Transmission Dynamics

This paper is concern with modeling cholera epidemic. Despite the advances made in understanding this disease and its treatment, cholera continues to be a major public health problem in many countries. Deterministic and stochastic models emerged in modeling of cholera epidemic, in order to understand the mechanism by which cholera disease spread, conditions for cholera disease to have minor and large outbreaks. We formulate a continuous time Markov chain model for cholera epidemic transmission from the deterministic model. The basic reproduction number (R0) and the extinction thresholds of corresponding cholera continuous time Markov chain model are derived under certain assumptions. We find that, the probability of extinction (no outbreak) is 1 if R0 < 1, but less than 1 if R0 > 1. We also carry out numerical simulations using Gillespie algorithm and Runge–Kutta method to generate the sample path of cholera continuous time Markov chain model and the solution of ordinary differential equation respectively. The results show that the sample path of continuous time Markov chain model fluctuates within the solution of the ordinary differential equation.


Introduction
Cholera is a virulent disease that affects both adults and children, and can kill within few hours if remain untreated.It is a severe water borne infectious disease caused by the vibrio cholerae bacterium (World Health Organization, 2010).Cholera bacterium is found in water or food sources that have been contaminated by feaces from a person who is infected with cholera.Cholera is most likely to be found and spread in places with inadequate clean water, poor sanitation, and inadequate hygiene.Therefore, one can get cholera by drinking water or eating food that are contaminated with the cholera bacterium (Ali et al, 2012).
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 worldwide.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 (World Health Organization, 2016).For example, in Tanzania the most recent information on cholera suggest that between August 2015 to April 2016, cholera resulted in over 24,000 cases and caused 378 deaths (World Health Organization, 2016).
There are several measures which have been suggested by (World Health Organization, 2010) to prevent cholera.These measures are environmental sanitation, provision of clean water, provision of education on the effect of cholera (World Health Organization, 2006).Despite various efforts made against cholera, the disease still persist in the community and is considered as a very threatening disease in many countries.
A number of studies have been done to study the spread of infectious diseases in deterministic and stochastic models setting (Dimi et al, 2016, Edward & Nyerere, 2015, Einsenberg et al, 2013, Khan et al, 2015, King et al, 2008, Longini et al, 2007, Neilan et al, 2010, Njagarah & Nyabadza, 2015, Panja & Mondal, 2014, Tuite et al, 2013, Marwa et al, 2018).Some of the studies consider compartment I for infected individuals as a single compartment instead of splitting it into two heterogeneous groups that is symptomatic and asymptomatic infected individuals to study the transmission dynamics of cholera epidemic.Also, continuous time Markov chain models, were ignored in their study.
Continuous time Markov chain models help to determine the probability of an outbreak or disease extinction through the use of branching process theory.In continuous time Markov chain we can study the stochastic thresholds for the disease extinction or major outbreak (Allen & Lahodny, 2012).The stochastic thresholds in continuous time Markov chain models are similar to the basic reproduction number, but these stochastic thresholds usually depend on the initial number of infectious individuals and the initial level of the environmental contamination.However, the stochastic thresholds are usually studied from the theory of multitype branching process.Berriman et al. (2013) developed a multi-group semistochastic model to describe Salmonella dynamics.In their study did not consider the probability of a major outbreak or host demography.The continuous time Markov chain for cholera epidemic developed in (Lahodny et al, 2015).The model accounted for host demography, direct and indirect transmission, replication of free-living pathogens in the environment and removal of free-living pathogens by natural death or environmental decontamination.The limitations to their model, the infected individuals were termed as a homogeneous group and also, they did not account for removal of free-living pathogens from the environment due to ingestion by hosts.
In this paper we, thus extend the deterministic model developed in (Marwa et al, 2018) by formulating an equivalent cholera continuous time Markov chain model.The formulated cholera continuous time Markov chain (CTMC) model will be analyzed using branching process approximation and probability of extinction.Also, numerical simulation will be carried using Gillespie algorithm.The advantage of using continuous time Markov chain (CTMC) model over deterministic model in epidemic modelling is that the CTMC model yields an estimate of the probability of extinction and outbreaks of cholera disease (Zevika & Soewono, 2018).Also, the measured trajectories behave as predicted.Therefore, the overall purpose of this paper is to study the transmission dynamics of cholera using continuous time Markov chain.
In the next section, we present a deterministic cholera epidemic model with water treatment, which is a model formulated by (Marwa et al, 2018).In section 3, the corresponding cholera continuous time Markov chain model is formulated basing on the assumptions of the ODE model.In section 4, the multitype branching process is discussed.In section 5, the numerical results from the deterministic and stochastic simulations are presented, discussed and compared.Finally, the last section concludes our paper  (Marwa et al, 2018) 1.From the description of the dynamics of cholera as shown in Figure 1, we have the following system of nonlinear ordinary differential equations (Marwa et al, 2018).

Reproduction Number, R 0
The key parameter in epidemiology is the basic reproduction number R 0 , which is defined as the average number of secondary infectious cases transmitted by a single primary infectious cases introduced into a whole susceptible population (Cohn et al, 1991).This parameter is useful because it helps determine whether an infectious disease will spread within the population or not.To compute R 0 , the next generation matrix approach is used as in (Van den Driessche & Watmough, 2002).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 (t), I a (t) and B(t) in the model (2.1).Hence, the basic reproduction number obtained in (Marwa et al, 2018) was as follows: where D = δκµd + δκµ 2 + δκµr 1 + δκdr 2 + δκµr 2 + δκr 1 r 2 + δκµϕ + ϕκµ 2 + κµϕr 1 + κϕdr 2 + κµϕr 2 + κϕr 1 r 2 .

Cholera Continuous Time Markov Chain Model
From the deterministic model (2.1), a cholera continuous time Markov chain (CTMC) model is formulated.The cholera continuous time Markov chain is a time-homogeneous process with the Markov property.By considering the spread of epidemic disease, clearly it is more realistic to keep track of where the epidemic disease is at minor, large outbreak or endemic at any continuous time t ≥ 0 as opposed to only where the epidemic disease is after n steps under various conditions.
The cholera CTMC model gives an indication to the stochastic variability inherent in birth, death, and transmission which are tracked in the deterministic model.However, cholera CTMC model provide enough information for the given population dynamics when the population sizes are small.Also, these birth, death, transmission and other transitions are expressed in terms of probabilities.
The cholera CTMC epidemic processes are defined on a continuous time scale, t ∈ [0, ∞), but the states that denotes the concentration of vibrio cholerae in the environment at time t.However, B(t) is not a compartment occupancy as the other one are so its transition is not considered in the formulation of transitions.Now suppose that, the population is homogeneously mixed so that each susceptible individual has equal chance of acquiring cholera through consuming water which are contaminated with vibrio cholerae in the reservoir or environment B(t).The susceptible individuals change their classification and belong to class I a (t) or I s (t).Class I a (t) and I s (t) change to R(t), where I a (t), I s (t) will remain indefinitely, independent of the past.The length of time spent in state I a (t), I s (t) or R(t) is continuously and strictly positive.
In order to develop the transition probability matrix p s (u) of possible states for cholera CTMC model, we denote the change for each variable during time △t as The lower case (s, i s , i a , r) denotes the values of the discrete random variables from the set {0, 1, 2, 3, • • • , N}.The transition probabilities together with stochastic process are defined for a small time interval of length △t > 0 as p s ((s, i s , i a , r), (s + e 1 , i s + e 2 , i a + e 3 , r + e 4 )) (△t) = (s + e 1 , i s + e 2 , i a + e 3 , r + e 4 )|(S (t), I s (t), I a (t), R(t)) = (s, i s , i a , r)).The transition probabilities depend on the time between events △t but not on the specific time t, known as a time-homogeneous process.However, when given the current state of the process at time t, the future state of the process at time t +△t, for any △t > 0, does not depend on time prior to t, this is termed as Markov property (Allen, 2017).Table 2, shows the changes, , associated with four events, birth, death, infection and recovery.
rate of recovery of I a Basing on the comparison purposes, the transition probabilities are defined in terms of the rates in the cholera model (2.1) as follows.

Kolmogorov Differential Equations
From Equation (2.5), we can derive differential equations for the transition probabilities.These differential equations can be referred to as the backward or forward Kolmogorov differential equations (Allen, 2017).The forward Kolmogorov differential equations, are also known as the Master equations, which are used to predict the future dynamics (Allen, 2017).Therefore, from Equation (2.5), we compute (p s (s, i s , i a , r) (t + △t) − p s (s, i s , i a , r) (t)) △t , and letting △t → 0, leads to the system of forward Kolmogorov differential equations: (2.6) It should be noted that both b and β may be constant or time dependent.Then by assuming that T j is the time of the j th jump, then the time of the ( j + 1) th jump is, where W j is the interevent time.Now, for x(T j ) = (e 1 , e 2 , e 3 , e 4 ) T , and n = e 1 + e 2 + e 3 + e 4 , we get the interevent time as where U is a uniformly distributed random variable on the interval [0,1] ( , is the mean for the exponentially distributed random variable W j .

Multitype Branching Process
The multitype branching process helps to approximate the dynamics of disease of nonlinear continuous time Markov chain model near the disease free equilibrium (dorman et al, 2004).In this paper, we use the idea of multitype branching process to determine the stochastic threshold.

Branching Process Approximations
In this paper, we also study the stochastic behaviour near the disease free equilibrium.This helps to know whether an epidemic lead to major outbreak or minor outbreak for cholera CTMC model when few infectious individuals are introduced into the population.This is possible by using branching process theory and other techniques from the probability generating functions.It should be noted that when states I s (t), I a (t) and B(t) reaches an absorbing states, then I s (t) = I a (t) = 0, and B(t) = 0.This implies that, the disease transmission stops (Allen, 2017).
However, the branching process is termed as the linear approximation of the SIR stochastic process when near the disease free equilibrium.Hence, when there is few initial infectious individuals, the branching process tend to hit zero or grows exponentially (Allen, 2017).The main two features which are captured in the branching process approximation for the continuous time Markov chain model near the disease free equilibrium are: 1.The major outbreak occurs when infectious individuals increases to a large number of cases and 2. The minor outbreak occurs when there is few additional cases above the initial number of cases.
Furthermore, it should be noted that the branching process is the best approximation for continuous time Markov chain model when the number of susceptible individuals are large (Allen, 2017).
The branching process is a birth and death process for I s (t), I a (t) and B(t), therefore the variable R(t) will not be considered in the branching process approximation and the susceptible individuals are assumed to be at the disease free state as follows Both infectious individuals I s (t) and I a (t) can produce secondary infectious hosts through shedding vibrio cholerae to the environment and hence, susceptible individuals will be at high risk of getting infection.Consider the cases; r 1 I s (t) be recovery rate of symptomatic infected individuals, dI s (t) death rate of symptomatic infected individuals due to disease, r 2 I a (t) recovery rate of asymptomatic infected individuals, and bN the birth rate or recruitment rate.The process start with few infectious individuals (i.e., symptomatic and asymptomatic infected individuals).
The following are the three important assumptions to underlie the branching process approximation (Allen, 2017).
1.Each infectious individual behaviour is independent from other infectious individual, 2. Each infectious individual has the same probability of transmitting an infection and same probability of recovery, 3. The susceptible population is very large.
From the assumptions with see that assumption ( 1) is feasible if a number of infectious individuals introduced into a large homogeneously mixed population is small.Where as assumption ( 2) is said to be feasible in a homogeneously mixed population with constant transmission and recovery rates, β, r 1 and r 2 (Allen, 2017).

Probability of Extinction
The probability of extinction can be studied by using two probability generating functions (pgfs).
1.The first probability generating function applies to each infectious individual known as the offspring pgf, and the second one 2. applies to the entire infectious class I s (t) and I a (t) at time t.
The term offspring, in this paper will help to describe susceptible individuals that become infectious due to indirect transmission, that is through the focal-oral route of contaminated food or water caused by poor sanitation (World Health Organization, 2010).The offspring pgf has the form: where p k is the probability of one individual generating k new individuals of the same type i.e., one infectious individual generate k infectious individuals and s is the extinction probability.
The probability generating function has some properties that are useful in the analysis, when f (1) = 1, then the mean number of offspring generated from one individual is defined as Referring to the nature of our model, the approximation of the continuous time Markov chain near the disease free equilibrium leads to the multi type branching process in the three dependent variables I s (t) , I a (t) and B(t).The three offspring pgfs are defined; two for infected individuals and one for vibrio cholerae concentration in the environment, where each offspring has the general pgf form: For each infected individual we have two events, that is infection and recovery.Therefore, the offspring pgf for I s (t) i.e., I s (0) = 1 and B(0) = 0 is given by From Equation (3.0), each symptomatic infectious individuals recovers with the probability r 1 r 1 + α 1 and each symptomatic infected host shed vibrio cholerae on the environment with the probability , where the term s 1 implies that each symptomatic infected individual generates an individual with vibrio cholerae (s 2 raised to the power of one) and they remains infectious (s 1 , s 2 raised to the power one).
For the case of asymptomatic infected I a (t), To find the probability of extinction (no outbreak), we compute the fixed points of the system (w 1 , w 2 ) ∈ [0, 1] of the offspring pgfs for the infectious classes, that is, we solve f 3 (s 3 , s 4 ) = s 3 , f 4 (s 3 , s 4 ) = s 4 , e.g., (Allen & Van den Driessche, 2013).We have the offspring pgf for I a (t) i.e., I a (0) = 1 and B(0) = 0. Hence In Equation (3.2), each asymptomatic infectious individuals recovers with the probability r 2 r 2 + α 2 and each asymptomatic infected shed vibrio cholerae on the environment with the probability , where the term s 3 implies that each asymptomatic infected individual generates an individual with vibrio cholerae (s 4 raised to the power of one) and they remains infectious (s 3 , s 4 raised to the power one).
Also infects individuals with the probability The two solutions for this system are (1, 1) and w 3 , w 4 , where ) , ) , and It should be noted that, the fixed point (w 1 , w 2 , w 3 , w 4 ) exists in the unit square provided that R 01 , R 02 > 1.Also, w 1 and w 3 have the biological interpretation, that is from the beginning of one infectious hosts, there is no outbreak if the infectious individuals recovers with the probabilities r 1 r 1 + α 1 and r 2 r 2 + α 2 or there is no successful transmission to susceptible cholera with probability ) and Then, the probability of no outbreak is 1, if R 01 < 1 and R 02 < 1 respectively.Also less than 1, if R 01 > 1 and R 02 > 1 respectively.Given that I s (0) = i s , I a (0) = i a , and B(0) = b.From the independent assumptions of branching process approximation, the probability of major and minor outbreak for symptomatic infected individuals are given by and where i s = I s .
Also, for asymptomatic infected individuals we have and The branching process theory when applied to infectious disease provide a clear method for the prediction of disease outbreaks.Whittle (1955) as one of the investigator of this theory showed that, the basic reproduction number R 0 , computed from the deterministic model yields the same threshold as the stochastic Markov chain SIR model, when given an initial infectious individuals i related to the population size N. Whittle (1955) also, stated that when R 0 < 1, it means no outbreak, while when R 0 > 1 there is an outbreak.Thus, the probability of a minor outbreak is ) i and the probability of major outbreak is 1 − ) i .These estimates for extinction are more accurate for a large susceptible population size N and when there is few infectious individuals i.Thus

Numerical simulation and Discussion
The model formulated in this paper involves multivariate processes, so it is difficult to find analytical solution for the transition probabilities from the forward Kolmogorov differential equation.Therefore, a numerical method for simulating the sample path of CTMC model formulated in this paper is to use Gillespie algorithm.This method was proposed by Gillespie (Gillespie & Marshall, 1977).Also, in order to numerically simulate the change in state, we need to have this uniform random numbers, u i ∈ U[0, 1] which are required for each change.Here u i stands for the particular interevent time and event.While i takes the rates for birth, death and movement between compartments.
Furthermore, to numerically compute the sample path of cholera continuous time Markov chain (CTMC) model, we need to employ a fact that interevent time has an exponential distribution.This is deduced from the Markov property, due to this an exponential distribution has a property called memoryless (Allen, 2008).The formulation of cholera CTMC model from the deterministic model (2.1) in this study is based on three rates: birth, death and movement between compartments.The birth, death and movement rates are defined as follows: the leaving compartment at a rate of λ m1 (k) and k → k + 1 is the entering compartment at a rate of λ m2 (k).The rates λ b (i), λ d (i), λ m1 (k) and λ m2 (k) denotes birth, death and movement of individuals respectively.The total rate λ is given by When an event occurs, it will be a birth with probability While an event dies with probability .6)The movement between individuals occur with probability Also, we apply a uniform random number on [0,1] for numerical computation of the interevent time.Now, let U be a uniform random variable on [0,1].Then the interevent time T i,k satisfies

Conclusion
The intent of this paper is to formulate a cholera continuous time Markov chain model from the deterministic model.Many papers of relevance to stochastic epidemic modelling omitted continuous time Markov chain.This paper provide enough information on cholera CTMC model and other topics in CTMC epidemic modelling.We apply the theory of multitype branching processes in order to estimate the probability of disease extinction or major outbreak.Also, we illustrate the effects of infectious individuals shedding vibrio cholerae into the environment and the environmental decontamination on the probability of disease extinction or major outbreak.From this study we find that, the probability of extinction (no outbreak) is 1, if R 0 < 1, but less than 1, if R 0 > 1 and when I s = I a = 0, and B = 0 the epidemic stops.This implies that an absorbing state has been reached and the process stops.However, numerical simulation is carried using Gillespie

Figure 1 .
Figure 1.The interaction of cholera epidemic transmission dynamics between compartments are shown .9) From the idea of Equations (3.4) to (3.9), four sample paths of the S I s I a R − B continuous time Markov chain model are graphed in Figures 2 to 5. Also, we numerically solve the ordinary differential equations (ODEs) model (2.1) using the fourth order Runge-Kutta method.A stable endemic equilibrium exists in the ordinary differential equations (ODEs) model if R 0 > 1 and the sample paths of S I s I a R − B continuous time Markov chain model fluctuates around this endemic equilibrium.The parameter values used in plotting the sample path of Figures 2 to 5 are described in Table 3, with the initial conditions S (0) = 997, I s (0) = 1, I a (0) = 1, R(0) = 1, N = 1000, B = 1, and △t = 0.01.The results show that the sample path of cholera continuous time Markov chain (CTMC) model fluctuates within the solution of ordinary differential equations as seen in Figures 2 to 5.

Figure 2 .
Figure 2. The sample path of the symptomatic infected individuals for CTMC model is graphed with the deterministic solution (dashed curve).The values are △t = 0.01, I s (0) = 1, N = 1000

Figure 4 .
Figure 4.The sample path of the susceptible individuals for CTMC model is graphed with the deterministic solution (dashed curve).The values are △t = 0.01, S (0) = 997, N = 1000 algorithm and Runge-Kutta method to simulate the sample path of S I s I a R − B continuous time Markov chain model and the solution of ordinary differential equation respectively.We find that, the sample path of S I s I a R − B continuous time Markov chain model fluctuates within the solution of the S I s I a R − B ordinary differential equation model.

Table 1 .
. In their model, S (t) denotes number susceptible individuals, I s (t) denotes the number of symptomatic infected individuals, I a (t) denotes Parameters of the cholera model