An Optimal Cost Effectiveness Study on Zimbabwe Cholera Seasonal Data from 2008–2011

Incidence of cholera outbreak is a serious issue in underdeveloped and developing countries. In Zimbabwe, after the massive outbreak in 2008–09, cholera cases and deaths are reported every year from some provinces. Substantial number of reported cholera cases in some provinces during and after the epidemic in 2008–09 indicates a plausible presence of seasonality in cholera incidence in those regions. We formulate a compartmental mathematical model with periodic slow-fast transmission rate to study such recurrent occurrences and fitted the model to cumulative cholera cases and deaths for different provinces of Zimbabwe from the beginning of cholera outbreak in 2008–09 to June 2011. Daily and weekly reported cholera incidence data were collected from Zimbabwe epidemiological bulletin, Zimbabwe Daily cholera updates and Office for the Coordination of Humanitarian Affairs Zimbabwe (OCHA, Zimbabwe). For each province, the basic reproduction number () in periodic environment is estimated. To the best of our knowledge, this is probably a pioneering attempt to estimate in periodic environment using real-life data set of cholera epidemic for Zimbabwe. Our estimates of agree with the previous estimate for some provinces but differ significantly for Bulawayo, Mashonaland West, Manicaland, Matabeleland South and Matabeleland North. Seasonal trend in cholera incidence is observed in Harare, Mashonaland West, Mashonaland East, Manicaland and Matabeleland South. Our result suggests that, slow transmission is a dominating factor for cholera transmission in most of these provinces. Our model projects cholera cases and cholera deaths during the end of the epidemic in 2008–09 to January 1, 2012. We also determine an optimal cost-effective control strategy among the four government undertaken interventions namely promoting hand-hygiene & clean water distribution, vaccination, treatment and sanitation for each province.


Introduction
Cholera is still a burning problem in underdeveloped and developing countries causing morbidity and mortality. In Zimbabwe, one of the most severe cholera outbreaks occurred in [2008][2009], that had been attributed as the worst African outbreaks in terms of its high case fatality rate (CFR) and shorttime extensive spread in some provinces. The outbreak, beginning in Chitungwiza, had duration from August 2008 to July 2009, ultimately ended with 98,592 reported cases and 4,288 reported deaths [1]. These massive outbreaks happened mainly due to Zimbabwe's poor health care system, shortage of good-quality food and clean drinking water [2]. An economic crisis within this period accelerated the deterioration of the country's infrastructure, including a breakdown of basic municipal services (such as sewage treatment and water supply in many areas) and medical facilities [3].
The provinces of Zimbabwe experienced a total of 2101 cholera cases over the period, 17th October, 2009 to 30th June, 2011 [4,5]. The substantial number of cholera cases in some provinces, e.g. Manicaland, Mashonaland West, Masvingo, Midlands, etc., both during and after the epidemic in 2008-09, indicate a plausible presence of seasonal forcing in cholera incidence in some of the provinces.
Well strategic deployment of cholera intervention/interventions in Zimbabwe may reduce future cases and deaths, although the projected effect of available cholera interventions is debatable [6]. A lot of suggestions have come out for preventing the cholera outbreak in those regions. Many regional and international organizations suggest providing clean water, hand-hygiene (Soap) promotion and construction & promotion of sanitary systems. Other groups are arguing for the vaccination program, although some experts suggest that the effect of vaccination will be modest [7]. Several professionals have also recommended usage of rehydration therapy for mild infections (v10% bodyweight loss) and usage of antibiotics (Erythromycin, Doxycycline and Ringer Lactate) for severe cases (w10% bodyweight loss) to reduce morbidity [8], cost of productive time loss due to illness, and bacterial shedding [9]. With proper treatment of cholera cases, the CFR should remain below 1% [10]. However, in terms of cost effectiveness, cholera vaccination is by far most costly intervention [11] with US$1,658 to US$8,274 yields one DALY (Disabilityadjusted life year) and gaining that same year through promoting hand-hygiene need US$3.35, making hand-hygiene the cheapest among cholera interventions [12]. Even though, promoting hand-hygiene heavily depends upon the availability of clean water. So an optimal balance among different types of interventions may significantly reduce the number of cholera cases and deaths at a minimal cost. Thus, a well-coordinated effort and an effective response to control an outbreak are the most important tasks.
To control future epidemics, a good understanding of cholera transmission dynamics is crucial and mathematical models can be utilized as a potential tool [6,[13][14][15]. Some earlier studies on cholera are based on the assumption of constant transmission rate between human and bacterial population over time [6,[13][14][15] but in food or waterborne infections, the role played by temporal forcing is more subtle and interesting. There is strong evidence that the multi annual dynamics of cholera are interlinked with long-term environmental factors [16][17][18].
To capture the presence of possible seasonal pattern within the data of reported cholera cases, we include the periodicity factor in our model. With these backdrops, we modified the model proposed by Hartley et al. [15] to include periodic slow-fast transmission and fitted to Zimbabwe's weekly cholera seasonal data starting from 2008-2009 epidemics to June 2011. Daily and weekly data were collected from Zimbabwe epidemiological bulletin [4], Zimbabwe Daily cholera updates [5] and Office for the Coordination of Humanitarian Affairs Zimbabwe [19]. The basic reproduction number (R 0 ) carries information about the persistence of a disease [20,21]. It is inversely proportional to the mean age of (first) infection; greater it is shorter the generation time, and the disease transmission will be more explosive [21,22]. Aforesaid data was used to estimate R 0 in periodic environment, for all provinces across the country. To the best of our knowledge, this is probably the pioneering attempt to estimate R 0 in periodic environment using the real data set of cholera epidemic. We perform a statistical test suggested by Roger [23], using weekly cholera incidence data from each province to justify the presence of seasonal trend. We also study the existence of any underlying pattern of temporal forcing in slow-fast transmission rate with the seasonality in cholera incidence, as observed in some provinces. We provide forecasts of cumulative cases and deaths from the end of epidemic in 2008-09 to 1 st Jan 2012 for different provinces in Zimbabwe and study the optimal intervention strategy/strategies by minimizing the cost of different cholera interventions.

Basic model structure
We modify the existing model [15] assuming temporal variations in two types of transmission rates (slow and fast). The existing model [15] assumes constant human population size (birth and death rates are equal), neglects the cholera-related death rate and assumes life time natural immunity to cholera if recovered. We have modified these assumptions in our model by incorporating variable human population size, cholera-related death rate and the effect of natural immunity loss to cholera (as it is now proven fact that natural immunity to cholera varies from less than one year to two years [24]). Our basic model is a system of five differential equations (see Equation 1 & 2) describing how individuals can move between different states of susceptibility or infection with cholera.
We categorize the total human populations at time t (denoted by N(t)), into susceptible (S(t)), infected (I(t)), and recovered (R(t)) classes. A constant recruitment rate (P H ) to human population, which is the product of human birth rate (m b ) and initial entire human population size (N(0)), is considered.
Individuals die naturally at a rate m d . All newly recruited individuals are assumed to be susceptible.
A Recent study showed that freshly shed V. cholerae from human intestines are short-lived and hyper-infectious in nature [25]. It out competes other V. cholerae grown in vitro, by as much as 700-fold for at least the first 5 to 18 hours in the environment [25]. After the hyper-infectious stage, V. cholerae organisms lose their competitive advantage and become low-infectious. This hyperinfectivity is a key factor to understand the explosive nature of human-to-human transmission in cholera outbreaks. Based on this fact, we classify the bacterial populations in two states, one hyperinfectious (B H ) state for fast transmission and this hyper-infectious bacteria decay to become low-infectious (B L ) state after some time, which causes slow environmental transmission.
Susceptible individuals gain infection by consuming water contaminated with the low-infectious bacteria (B L ) and the high-infectious bacteria (B H ) at rates l L~b L (t)B L K L zB L and ))g are the rates of ingesting low-infectious and high-infectious V. cholerae bacterium from the contaminated water, which are assumed to be time periodic with period 52 weeks. b L0 , b H0 denote the minimum transmission rate of low and high infectious V. cholerae respectively from the contaminated water, and d denotes the amplitude of seasonality. Infected individuals either have a natural death (at a rate m d ) or die due to extreme loss of fluid from their body during infection (at a rate m c ) or recover naturally from cholera infection (at a rate c). Recovered individuals are immune to reinfection, but this immunity wanes over time and eventually returns to the susceptible stage (at a rate v). During the period of infection; infected individuals excrete V. cholerae into water reservoirs around them (at a rate j). Since this bacterium is coming directly from infected human intestines, it is in hyper-infectious stage and decays, ultimately leads to low infectious stage (at a rate x). The natural death rate of low infectious bacterium is d L . Based on the above-mentioned assumptions, we construct the following system of non-linear differential equations: where, Model parameters and their interpretations with some parameters' base values, taken from previous studies, are given in Table S1. A flow diagram of the basic cholera model (1) is also given in Figure 1.

Source of data
Daily and weekly cholera incidence data for each province of Zimbabwe have been collected from Zimbabwe epidemiological bulletin [4], Zimbabwe Daily cholera updates [5] and Office for the Coordination of Humanitarian Affairs Zimbabwe [19] starting from 2008-2009 epidemics to June 2011. The data at the beginning of the epidemic are quite noisy. To smoothen out the initial fluctuations, the data is converted to weekly reported cases by aggregating daily case reports for the entire duration of the epidemic. Table 1 contains information about the starting points and the end points of the data for each province. Total numbers of data points vary across the provinces.

Model calibration
To calibrate the basic model (1), we have considered the weekly reported cholera cases and deaths from each province of Zimbabwe starting from 2008-2009 epidemics to June 2011. The model is fitted to the cumulative number of cases and deaths obtained from the weekly counts in each province.
The key parameters estimated from the data are the average transmission rate of hyper-infectious Bacterium (b H0 ), the average transmission rate of low-infectious Bacterium (b L0 ), the amplitude of seasonality (d), the mortality rate of human due to cholera infection (m c ) and the excretion rate of cholera infected individual (j). It is not realistic to assume the entire population of a province to be susceptible to cholera, as outbreaks generally occur in that part of the province where the basic amenities like proper drainage system, clean water and food are lacking. So, we first estimate the initial number of susceptible, infected and recovered human populations from the data by bounding the initial total human population size (N(0)) by the total population size of the province. The initial concentrations of hyper-infectious (B H (0)) and lowinfectious (B L (0)) bacterium are estimated from the data. Our work also involves estimation of initial reported cases (C(0)) and deaths (D(0)) from the data since in some provinces the exact reported cases and deaths from the beginning of the epidemic were unknown due to reporting delays.
The cumulative cases and cumulative deaths from the cholera model (1) are given by: whereh h[R d contains all the unknown variables of the model (1). We have n observations (cumulative cases and cumulative deaths) from our data at n different weeks t i as Y (t i )~(C(t i ),D(t i )) T , where i~1,2,:::,n and t i is the i th week in our data. We assume independent Gaussian prior specifications forh h: h j *N(n j ,g 2 j ), j~1,2,:::,d ð4Þ Let e be the error when fitting cumulative quantities (C(t,h h),D(t,h h)) from the model (1) to the observed data. Then e follows independent Gaussian distribution having unknown variance s 2 i.e. e*N(0,Is 2 ). For the error variance a Gamma distribution is used as a prior for its inverse: where the prior parameters S 2 0 and n 0 in (5) can be interpreted as the prior mean for s 2 and the prior accuracy as imaginary observations.
We construct the sum of squares function as: Posterior distribution of the model unknown variableh h is generated using Delayed Rejection Adaptive Metropolis algorithm (DRAM) [26,27] with an initial burn of 100000 iterations. MCMC toolbox in MATLAB written by Marko Laine [28] was used to estimate the unknown variableh h for the model (1). Geweke's Zscores [29] were examined to ensure the chain convergence.
The advantage of using the cumulative over the weekly number of new cases in model calibration is that the former smoothes out known reporting delays on weekends and national holidays [14,30].

Seasonality
To justify whether the existence of any kind of seasonal forcing influence the number of cholera incidence in Zimbabwe provinces, a suitable statistical testing procedure is very much needed. The weekly data from 18th August, 2008 to 30th June, 2011 is taken into account for this purpose. We follow the test procedure suggested by Roger [23].
The entire span of almost 3 years is divided into 52 classes, corresponding to the 52 weeks of an year (1 st week starting from 18 th August, 2008) and the total number of cholera cases in week i (i~1,2,::,52) is denoted by N i . The probability that any one event belongs to i th class is P i , where , i~1,2,::,52: where, q i denotes the frequency for class i under the null hypothesis, s i~s in (2pi=52), c i~c os (2pi=52), a and b are the parameters of the model (7). H 0 : a~b~0 indicates the absence of seasonality and H 1 : a=0 or b=0 indicate the seasonality in cholera incidence. The test statistics for testing H 0 is of the form Where, and n~X i N i : The test statistic R is asymptotically distributed as chi-square with 2 degrees of freedom.

Estimating reproductive numbers in periodic environments
The Model (1) has a unique disease free equilibrium given by: Following [31], we calculate the matrix of new infection from our system (1) as: Let, Y (t,s), t §s be the evolution operator of the linear $periodic system for all t §s and Y (s,s)~I, where I is the 3|3 identity matrix. Let C $ be the ordered Banach space of all $-periodic functions from R to R 3 which is equipped with maximum norm E:E ? and the positive cone , for all t in R}. Consider the following linear operator L : Following Wang and Zhao (2008) [31], we call L the next infection operator, and define the basic reproduction number (R 0 ) as: where r(L) is the spectral radius of the operator L defined in Equation (11). Motivated by the concept of the partial reproduction numbers defined by Mukandavire et.al. [14], we similarly define two partial reproduction numbers R l and R h in periodic environment. The subscripts l and h correspond to low infectious and high infectious transmission, respectively.
Basic reproduction number (R 0 ) is the sum of two partial reproductive numbers-the one is arising from the contact between the human and low-infectious bacteria, which we denote as R l and the other one arising from the contact between the human and hyper-infectious bacteria, denoted as R h . Using Lemma 1 given in Appendix S1 and the estimated parameter values (Table S3  and Table S4), we numerically estimate R 0 , R l and R h for each province. Details of the derivation procedure of R 0 , R l and R h are given in Appendix S1. For uncertainty, we draw 95% confidence interval around the estimated values. The following procedure is applied to derive the 95% confidence intervals for R 0 , R l and R h , respectively.
We draw a sample of size n (n~1000) from the posterior distribution ofh h (set of model variables, which are estimated) using simple random sampling without replacement (SRSWOR) scheme. The posterior distribution ofh h is depicted in Figure S1 and Figure S2. For each of the sample values ofh h, we estimate numerically the value of R 0 (using Lemma 1, Appendix S1). Thus, a vector of size n, is generated for R 0 . Curtailing the lower 2:5% and the upper 97:5% observations from the ordered vector of R 0 , we obtain the 95% confidence interval for R 0 . Applying the similar procedure we draw 95% confidence intervals for two partial reproductive numbers, R l and R h , respectively.

Projection of future cases and deaths
We project the number of cholera cases and deaths from the end of epidemic in 2008-09 to January 1, 2012. For uncertainty, we derive 95% credible intervals around the estimates of future projected cases and deaths. To predict the number of cases and deaths for a particular province, we simulate the cholera model (1), using the known & estimated parameters (Table S1 and Table  S3) and demographic parameters (Table S4), up to the end of epidemic in 2008-09, in that region. We obtain different demographic variables of human & pathogen (S, I, R, B H and B L ), new cases and deaths corresponding to the end of epidemic in 2008-09. Using this information from the previous simulation as initial conditions and parameter values from Table S1 and  Table S3, we simulate the model (1) to obtain predicted cases and deaths from the end of the epidemic in 2008-09 to January 1, 2012.
We used the following procedure to derive the 95% confidence intervals for projected cases and deaths. For each of the sample value ofh h (see, section-Estimating reproductive numbers in periodic environments), we predict the number of cases and deaths using the above procedure. Thus, two vectors, each of size n, are generated for predicted cases and deaths, respectively. Curtailing the lower 2:5% and the upper 97:5% observations from the ordered vector of predicted cases and deaths we obtain the 95% confidence intervals for projected cases and deaths, respectively.

Model with different cholera interventions
Effect of four different types of cholera interventions namely hand-hygiene promotion & clean water supply, treatment using oral rehydration therapy & antibiotics, vaccination and sanitation are studied. We assume that hand-hygiene (soap) & clean water will reduce bacterial ingestion by a fraction (1{h(t)), where h(t) is the relative rate of reduction in bacterial ingestion per week using hand-hygiene & clean water supply. Vaccinated population is increased by a proportion p(t)s of the susceptible individuals, who are successfully vaccinated, where p(t) is the per week vaccination rate and s is the vaccine efficiency. Vaccinated population is decreased due to the waning of vaccine based immunity (at a rate e) to become susceptible again and die (natural deaths) at a rate m d . We assume that a proportion a(t) of the infected individuals receive treatment by oral rehydration salt (for v 10% body weight loss) and by antibiotic (for w10% body weight loss) per week. Natural recovery rate of the treated person increases by relative rate of recovery l. Since, excretion can be affected by the use of antibiotics [32,33], the relative rate of shedding is reduced by a fraction y among the proportion of the infected individuals who receive antibiotic at a rate a(t) per week. Now proper sanitary and drainage system will prevent the human waste to contaminate the nearby water reservoirs. Invariably, sanitation will reduce excretion rate of human that contaminate the nearby reservoirs by a fraction (1{s(t)), where s(t) is the rate of reduction in human shedding per week by construction and promotion of sanitation. In Zimbabwe, 80% of the total populations have access to improved water source and 40% of the total populations have access to proper sanitation facilities [1]. Therefore, we assume maximum percentage reduction in bacterial ingestion rate through hand hygiene & clean water supply and reduction in human shedding possible by promoting sanitation in a week to be 80% and 40%, respectively. We also assume that maximum 70% of the infected individuals receive proper treatment in a week and maximum vaccination coverage possible in a week is about 35% of the total susceptible population [34]. Effects of hand-hygiene & clean water, vaccination, treatment, sanitation and their different combinations are projected from the end of epidemic in 2008-09 to January 1, 2012.
System of non-linear differential equations representing the effect of different interventions on our basic model (1) is given as follows: Intervention parameters and their interpretations with some parameters' base values taken from earlier studies are given in Table S2. A flow diagram of the intervention model is depicted in Figure 2.

An optimal intervention strategy
To determine the optimal intervention strategy/strategies (which reduce the number of cases and deaths projected from the end of the epidemic in 2008-09 to January 1, 2012, at a minimal cost), we define the following cost function: where interventions are applied for T weeks. First term in the right-hand side of (14) represents the cost of cholera-related deaths and remaining terms are costs associated with the implementation of different interventions. Nonlinear terms in the objective function J represent the costs of interventions in emergency situations. A, B, C, D, E, F , G, H and K are fixed cost coefficients, given in the Table 2.  Table 2) respectively. Y intervention and Y cost year are in the same base 1982~100, data of annual consumer price index (US) were collected from U.S. Bureau of Labor Statistics [35].
Our goal is to minimize the objective function J with respect to different control parameters h(t), p(t), a(t) and s(t) to determine an optimal intervention combination. This is a dynamic control problem and is solved directly by using the Pontryagin's Maximum Principle [36] and the method of steepest decent [37]. Minimization procedure of the objective function J is briefly described in Appendix S1.
The average coverage percentages (per week) of hand-hygiene (soap) & clean water distributions, treatment and sanitation are estimated using the following formulas:ĥ where,ĥ h,â a andŝ s, denote the average coverage percentages (per week) of hand-hygiene (soap) & clean water distributions, treatment and sanitation, respectively. h c (t), a c (t) and s c (t) are the optimal rates of hand-hygiene (soap) & clean water distribu-  tions, treatment and sanitation, respectively, for which the cost function J (see equation (14)) is minimum. T denotes the total number of weeks during which an intervention is applied. Total coverage percentage of vaccination is estimated using the following formula: Total vaccination coverage~V where, V (0) is the initial number of vaccinated individuals and S(t) is the number of susceptible individuals at week t. s is the vaccine efficiency (Table S2). p c (t) is optimal vaccination rate which minimizes the cost function J (see equation (14).
Cost per averted case for an intervention is calculated using the following formula:

Total cost of an intervention Total number of cases averted by that intervention
A 95% confidence interval for each of the following quantities are obtained following the same technique as explained in sections Projection of future cases and deaths and Estimating reproductive numbers in periodic environments: (1) cases that occurred in spite of applying an intervention, (2) total cost of an intervention and (3) cost per averted case.

Results
Cholera model fitting for the cumulative reported cholera cases and deaths are depicted in Figure 3 and Figure 4, respectively. A comparison between weekly reported cholera cases & deaths from each province with the model solution are shown in Figure 5(A, B, and C) and Figure 6(A, B, and C), respectively. The estimated model parameters, including human and pathogen demographic parameters, for each province are given in the format [estimate (95% CI)], (Table S3 and  Table S4). Plots for the posterior distributions of the estimated unknown variables of the cholera model (1) are given in Figure  S1 and Figure S2.
The contributions of low-infectious (R l ) versus high-infectious (R h ) transmission to R 0 vary widely. Our estimated values of R 0 , R l and R h are in good agreement with the previous estimates given by Mukandavire et. al. [14], for the provinces Harare, Mashonaland East, Mashonaland Central, Midlands and Masvingo but significantly differ in Bulawayo, Mashonaland West, Manicaland, Matabeleland South and Matabeleland North. In Mashonaland West and Manicaland, contribution of low-infectious (R l ) is higher than hyper-infectious (R h ) transmission to R 0 (see Table 3). Opposite trend is observed in the estimates of R l and R h given by Mukandavire et. al. [14] in these two provinces. The observed data points (available at some discrete time points over a time period, which varies across the study regions) are shown by blue circles while the solid lines depict the model solutions. The cumulative cholera cases from the model are plotted for each day of the time period (from the start to end week for the observed cholera data) using parameter values and initial conditions from Table S3 and Table S4.   Table S3 and Table S4.   Table S3 and Table S4.   Table S3 and Table S4. doi:10.1371/journal.pone.0081231.g006 Estimates of R 0 in Bulawayo (0.1022), Matabeleland South (0.7914) and Matabeleland North (0.0541) are all found to be below the unity, which are drastically different from the estimates given by Mukandavire et. al. [14].
To justify the existence of any seasonal trends in cholera incidence data of different provinces, the values of the test statistic R, defined in (8), are calculated using weekly data from 18th August, 2008 to 30th June, 2011. The values of R and corresponding p-values for each province are given in Table 4. Significant seasonal trends in cholera incidence data are observed in Harare, Mashonaland West, Mashonaland East, Manicaland and Matabeleland South. Among these five provinces, Harare and Manicaland exhibit highly significant seasonal trend, as confirmed by the corresponding p-value of the test (v 0.01).
Our basic model (1) without any interventions, projects 6340(5565{7264) cases and 271(238{309) deaths due to cholera in Zimbabwe between the end of epidemic in 2008-09 to January 1, 2012. Table 5 and Table 6 contain the predicted total cholera cases and deaths for the provinces during that time interval. Among the ten provinces, in Mashonaland West and Mashonaland Central the most numbers of cholera cases and deaths are predicted. In Mashonaland West, total 3118(2788{3463) cases, and 135(120{151) deaths are predicted during the mentioned period. This is about 49% of the total predicted cases and about 50% of the total predicted deaths in all provinces of Zimbabwe. In Mashonaland Central, total 987(858{1161) cases, and 33(29{39) deaths are predicted during the mentioned period. This is about 16% of the total predicted cases and about 12% of the total predicted deaths in all provinces of Zimbabwe.
To justify the predictive performance of our basic model (1), we compare the predicted cumulative cases and deaths from the end of the epidemic in 2008-2009 to January 1, 2012 with the reported cases and deaths' figures. The reported cumulative cases and deaths during the aforesaid time period are 2225 and 72, respectively [38], which is about 35%(31%{40%) of the model projected cases and about 27%(23%{30%) of the model projected deaths. Significant difference in the actual and predicted case and death figures may be attributed to the higher percentage of underreporting of cholera cases and deaths [39], that was not considered while making these predictions. According to WHO, the officially reported cholera cases represent only 5{10% of the actual number of cases those are occurring annually worldwide [39].
We found that, in the African region the countries report cholera cases more consistently than the other countries under WHO [39]. Also Zimbabwe's Integrated Diseases Surveillance & Response Technical guidelines list Cholera among the diseases that must be reported on a daily basis during epidemics to prevent avoidable illness and death [40]. Thus we may expect that the percentage of reported cases is higher in Zimbabwe than the worldwide statistics of under-reporting, although, we do not have the specific figures/numbers from literature depicting the actual percentages of under-reporting in Zimbabwe during the end of epidemic in 2008-09 to January, 1, 2012. The actual reported cases during that period are about 31%-40% of our model predicted cases. Hence this percentage may be considered as an estimate of reporting of cholera cases in Zimbabwe, which is much greater than the worldwide statistics (5%-10%).
It is already mentioned that Mashonaland West and Mashonaland Central are high-risk provinces in terms of cholera incidence between the end of epidemic in 2008-09 and January 1, 2012. Therefore, we discuss the results of different interventions and their layered combinations for these two provinces only.

Discussion
Our analysis suggests that, the routes of cholera transmission vary from province to province, that agrees with the findings of Mukandavire et. al. [14]. This heterogeneity in transmission dynamics may be due to the diverse geographic and climatic conditions across the country. A similar pattern in transmission dynamics is observed in Harare, Mashonaland West, Manicaland and Matabeleland South, where seasonality in cholera incidence was observed. In these provinces, slow transmission route is a dominating factor (b L0 w b H0 ) for cholera transmission (Table  S3). Earlier studies on cholera suggest that slow transmission route is more correlated to the climatic and environmental factors [17,18,[41][42][43][44][45][46][47] and is the main cause for seasonal dynamics of cholera [41,[48][49][50][51]. Unfortunately, due to lack of climatic data of Zimbabwe, we are unable to draw any quantitative inference, for example, whether inter-annual climatic variation in different provinces affects the transmission dynamics of cholera or not.
Estimate of human shedding rate (j) ( Table S3) in Bulawayo province differs from other nine provinces (one order of magnitude higher than for the other provinces). A possible reason for such difference may be due to the fact that the city receives portable water supply from five surface dams [52] and constantly suffers   from improper waste management system [53]. In spite of this; the national water supply agency of Zimbabwe (ZINWA) is not in charge in supplying water in Bulawayo [54].
Mukandavire et. al. [14] estimated R 0 for 2008-2009 cholera epidemics in Zimbabwe with constant transmission rate but to our knowledge, this is the first time that the basic reproduction numbers with periodic transmission rate are estimated for cholera epidemics in Zimbabwe or any other country. It is also to be noted that our data set contains much longer time scale (from the beginning of cholera epidemic in 2008-2009 to June 2011) than the earlier analyzed data [14]. It is already pointed out that the estimated values of R 0 in Bulawayo, Mashonaland West, Manicaland, Matabeleland South and Matabeleland North by Mukandavire et. al. [14] differ with our estimates. The disease dynamics of cholera may not be captured properly by assuming constant contact rate between human and bacterial populations over time, as it also depends on temporal forcing. Thus, the model with periodic environment is more appropriate than the previous studies. We believe that the prediction, thereby proposed, will be helpful for policy makers. Optimal cost effective study in Zimbabwe, from the end of epidemic in 2008-09 to January 1, 2012, suggests that, as a single intervention hand-hygiene & clean water supply is the most costeffective way to control a future cholera outbreak in those regions where slow transmission is the dominating factor for cholera transmission (see , Tables 8 and 9). In terms of cases and deaths reduction during epidemic hand-hygiene & clean water supply is by far the most effective individual intervention among the other single interventions (see Tables 5 and 6). This result is in good agreement with the observations by Andrews and Basu [6], where they argued that hand-hygiene & clean water distributions will avert more cases and deaths than treatment and vaccination during the epidemic in Haiti. Treatment is the most cost-effective in those regions where hyper-infectious transmission is the main factor for cholera transmission (see , Tables 8 and 9). This result is in well agreement with the previous observation of Naficy et. al. [55] on the control of cholera in sub-Saharan refugee settings. Treatment and hand-hygiene & clean water supply with any other intervention combination will also be cost-effective, which could avert thousands of cases and hundreds of deaths in Zimbabwe at a minimal cost. A synchronized, timely and efficient intervention will effectively reduce the severity of the disease and number of deaths. Our mathematical model and its prediction will help the public health authority of Zimbabwe for making suitable intervention strategies. We also believe that, a similar method can be applied for endemic and epidemic cholera outbreaks in other regions/countries as well, in particular, with seasonal patterns in disease transmission.