-
Notifications
You must be signed in to change notification settings - Fork 0
/
cholera-rainfall.tex
254 lines (211 loc) · 50.4 KB
/
cholera-rainfall.tex
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
\begin{fullwidth}
\chapter[Rainfall as a driver of epidemic cholera: comparative assessments of the effect of intra-seasonal precipitation events]{Rainfall as a driver of epidemic cholera:\\comparative assessments of the effect of\\intra-seasonal precipitation events}
\chaptermark{Model comparison assessment of rainfall as a driver of epidemic cholera}
\label{ch:cholera-rainfall}
The correlation between cholera epidemics and climatic drivers, in particular seasonal tropical rainfall, has been studied in a variety of contexts. Several mechanistic models have included rainfall as a driver of cholera transmission by focusing on two possible pathways: either by increasing exposure to contaminated water (\eg due to worsening sanitary conditions during water excess) or by increased water contamination with freshly excreted bacteria (\eg due to washout of open-air defecation sites or overflows). In this chapter, the explanatory power of these different model structures is assessed by formal model comparison using deterministic and stochastic models of the type susceptible-infected-recovered-bacteria (SIRB). The incorporation of rainfall effects is generalized using a nonlinear function that can increase or decrease the relative importance of large precipitation events. The modeling framework is applied to the daily epidemiological data collected during the 2015 cholera outbreak in Juba, South Sudan. This epidemic is characterized by a particular intra-seasonal double peak of incidence in apparent relation with particularly strong rainfall events. Results suggest that rainfall-based models in both their deterministic and stochastic formulations outperform models that do not account for rainfall. In fact, classical SIRB models are not able to reproduce the second peak, suggesting it was rainfall-driven. Moreover stronger support is found for rainfall acting on increased exposure rather than on exacerbated water contamination. Although these results are context-specific, they stress the importance of a systematic and comprehensive appraisal of transmission pathways and their environmental forcings when embarking on the modeling of epidemic cholera.
This section is adapted from:
\longfullcite{Lemaitre:RainfallDriverEpidemic:2019}, referred in the following as the postprint, and its supplementary information as \textsc{si}.%\footnote[][-4\baselineskip]{JL, JP-S and DP developed the numerical tools required for the comparison. JL and DP conducted the numerical analyses for the deterministic model while JP-S conducted the numerical analyses for the stochastic model. AR and DP designed the framework for this study. JFW provided the epidemiological data. CS analyzed the data and prepared the model input. All authors contributed to the discussion of the results and writing of the manuscript.}.
\end{fullwidth}
\section{Rainfall and the transmission of cholera}\label{sec:rainfall-cholera-transmission}
Two main exposure pathways fuel cholera transmission across endemic and epidemic settings.
First, an \textit{indirect} exposure occurs from consumption of water contaminated by untreated sewage%\footnote{discovered by John Snow during the 1854 Broad Street cholera outbreak\fullcite{Snow:ModeCommunicationCholera:1855}}
; rainfall and the ensuing hydrologic transport processes might play a role in water contamination, for instance through the washout of open-air defecation sites and sewage circulation in the environment.
Alternatively, \textit{direct}, or human-to-human exposure occurs when the bacteria is transmitted from an infectious individual directly to a susceptible one, for example via contaminated food or fomite. In this case, environmental factors do not play a major role, except for transmission changes due to behavior modifications.
The mechanism behind the transmission of cholera has been postulated to be a combination of environmentally-mediated and direct exposures\shortcite{Sugimoto:HouseholdTransmissionVibrio:2014,Bi:MicroscaleSpatialClustering:2016,Lessler:MeasuringSpatialDependence:2016,Rinaldo:ModelingKeyDrivers:2017}.
% Cholera and climate (STOP)
\paragraph{Cholera and Rainfall} The relationship between different climatic and environmental factors and the transmission of cholera has long been studied. Works linking cholera outbreaks to anomalies in the El Niño Southern Oscillation\shortcite[-2\baselineskip]{Colwell:GlobalClimateInfectious:1996, Pascual:CholeraDynamicsNinoSouthern:2000, Hashizume:DifferentialEffectIndian:2013} have paved the way for a new field in epidemiological research.
\begin{marginfigure}[-1\baselineskip]
\centering
\includegraphics{fig/cholera-rainfall.png}
\margincaption[Similarities between daily cholera incidence and rainfall in Haiti]{\footnotesize Daily cholera cases (red) and daily rainfall (blue) in Haiti from September 15, 2010 to October
16, 2011. It highlights the visual correlation between heavy rainfall event and cholera outbreaks. Adapted from \fullcite{Gaudart:SpatioTemporalDynamicsCholera:2013}.}
\label{fig:rain}
\end{marginfigure}
Previous studies have highlighted the role of climatic drivers on cholera dynamics, mostly focusing on climate change effects on disease spread\shortcite{Rinaldo:ModelingKeyDrivers:2017,Hashizume:EffectRainfallIncidence:2008,Magny:CholeraOutbreakSenegal:2012,Rodo:ClimateChangeInfectious:2013,Ramirez:NinoClimateCholera:2016,Vezzulli:ClimateInfluenceVibrio:2016} or on the impacts of spatial and temporal heterogeneities\shortcite{Reiner:HighlyLocalizedSensitivity:2012, Baker-Austin:EmergingVibrioRisk:2013, Vezzulli:OceanWarmingSpread:2013, Cash:CholeraShigellosisDifferent:2014,Escobar:GlobalMapSuitability:2015, Vezzulli:EffectsGlobalWarming:2015,Perez-Saez:ClimatedrivenEndemicCholera:2017}. The effect of temperature on cholera transmission has been described (though mainly through the lense of \textit{Vibros} survival and ecology in natural environments), but for rainfall, it remains to be fully elucidated, possibly due to the multiple ways in which it can influence transmission at the local and regional scales\shortcite{Rinaldo:Reassessment20102011:2012,Eisenberg:ExaminingRainfallCholera:2013,Baracchini:SeasonalityCholeraDynamics:2017}. Indeed, intense rainfall events have been shown to alter infection risk through a variety of potential mechanisms, including: flooding, leading to sewage contamination of water sources\shortcite{Ruiz-Moreno:CholeraSeasonalityMadras:2007, Hashizume:EffectRainfallIncidence:2008}; increased hydrologic transport-driven iron availability in environmental waters that enhances pathogen survival and the expression of toxins\shortcite{Lipp:EffectsGlobalClimate:2002,Faruque:SeasonalEpidemicsCholera:2005, Hill:ToxigenicVibrioCholerae:2011}; dry spells inducing persistent low water levels leading to increased use of unsafe water sources\shortcite{Rebaudet:EnvironmentalDeterminantsCholera:2013}; and crowding during strong flood events\shortcite{Reiner:HighlyLocalizedSensitivity:2012}.
Most countries where associations between rainfall and cholera risk have been studied experience endemic cholera transmission. Empirical studies have shown a range of correlations, both positive and negative, endowed with time lags ranging from weeks to months\shortcite{Ruiz-Moreno:CholeraSeasonalityMadras:2007,Emch:SeasonalityCholera1974:2008,Magny:CholeraOutbreakSenegal:2012}. In general, rainfall has been found to increase cholera transmission, but there is evidence of the inverse, possibly due to pathogen dilution\shortcite{Ruiz-Moreno:CholeraSeasonalityMadras:2007}. Such variability reflects the variety of potential mechanisms whereby rainfall may alter infection risk. Similarly, a clear empirical correlation between intense rainfall and enhanced transmission is found in several regions hit by cholera epidemics\shortcite[][but Haiti is treated specifically in the next chapter]{Magny:EnvironmentalSignaturesAssociated:2008,Magny:CholeraOutbreakSenegal:2012,Rebaudet:EnvironmentalDeterminantsCholera:2013,Rebaudet:CholeraCoastalAfrica:2013,Jutla:WaterMarkerMonitored:2013} (fig. \ref{fig:rain}). Results therein showed that at all spatial scales and locations examined, tropical storms were correlated with increased cholera incidence with lags of the order of a few days. As a consequence, accounting for the related forcing of dynamic models resulted in improved fits of reported incidence.
Properly incorporating the effects of rainfall in mathematical models of cholera transmission is thus paramount to discriminate among the above-mentioned alternative transmission pathways, thus unlocking a predictive framework to evaluate the potentially rainfall-sensitive efficacy of available intervention strategies in endemic and epidemic settings. This becomes critically important when evaluating the number of averted infections by deploying vaccines, as was done in the aftermath of the passage of Hurricane Matthew\shortcite{Pasetto:RealtimeProjectionsCholera:2017}, or considering optimal deployment in space and time.
The effect of rainfall on indirect cholera transmission has been accounted for in two main fashions in recent mathematical models. On one side, a \textit{contamination}-centered approach suggesting that bursts of infections could be linked to increased contamination of the water compartment\shortcite{Rinaldo:Reassessment20102011:2012}, as in the ECHO model presented in \textsc{Chapter~2}. This process conceptualizes the washout of open-air defecation sites by hydrologic transport. The same `transport' effect may be realized by sewer collectors' overflows. In fact, both mechanisms have the net effect of charging progressively the bacterial concentration in the water reservoir\shortcite{Codeco:EndemicEpidemicDynamics:2001}. Pathogens' loads are washed out from a hydrologic catchment enclosing human settlements and their infective individuals shedding bacteria. Therein, pathogen survival and thus the toxicity of their loads depend on hydrologic residence time distributions\shortcite{Rinaldo:Reassessment20102011:2012,Rinaldo:ModelingKeyDrivers:2017}. Such loads increase as a function of rainfall, which acts as a proxy of runoff volumes. The second approach is \textit{exposure}-centered and employs a rainfall-dependent exposure rate subsuming both pathogen availability and the probability of the ingestion of contaminated water during wet spells\shortcite{Eisenberg:ExaminingRainfallCholera:2013}. Although both approaches are physically plausible, they have not been compared directly on the same datasets within a formal statistical framework, which would allow to highlight their respective merits and further recommendations for their use in different settings.
In this chapter, the explanatory power of these different types of rainfall-driven mechanistic models applied to a cholera outbreak in South Sudan is compared. The link between rainfall and cholera during the outbreak recorded in Juba in 2015 when an intra-seasonal peak of cholera cases was recorded possibly in correspondence to intense precipitation events is quantitatively examined. The analysis of the lagged relationship between rainfall rates and revamped cholera incidence is addressed via dynamical compartmental models considered both in deterministic and stochastic versions incorporating both direct (human-to-human) and indirect (water-to-human) disease transmission, and rainfall effects on both contamination and exposure.
\section{Case study: the 2015 cholera outbreak in Juba, South Sudan}\label{sec:data sets}
\begin{figure}\centering
\includegraphics{fig_cholera-rainfall/Lemaitre_ACTROP_2018_42_R1_fig2.png}
\caption[Cholera cases and precipitation during the 2015 cholera epidemic in Juba][2\baselineskip]{Reported cholera cases (dots) and precipitation (bars) during the 2015 cholera epidemic in Juba. The timing of the vaccination campaign is highlighted in yellow.}\label{fig:report}
\end{figure}
In the past years, South Sudan had been struck by several cholera outbreaks\footnote{More details on the history of cholera in South Sudan are found in: \fullcite{Sciarra:MathematicalModelingCholera:2018}}. Here the analysis focuses on the 2015 outbreak in Juba, when a particular double peak of cholera cases was associated with a strong intra-seasonal precipitation event (fig.~\ref{fig:report}). Access was granted to epidemiological records for the 2014 and 2015 cholera epidemics that include daily cholera incident cases and hospitalization at the second-lowest administrative level (Payams), reported by the Ministry of Health of South Sudan. The cases in the 7 Payams that constitute the administrative area of Juba have been aggregated to obtain the reported time series for the county level. The population data for Juba county is taken from the South Sudan National Bureau of Statistics\cite[-4\baselineskip]{SSNBS:PopulationProjectionsSouth:2015}. Daily rainfall estimates for years 2014 and 2015 were obtained from the Climate Data Library\cite[-3.5\baselineskip][with a resolution of {0.1}$^\circ$, rainfall was averaged over the study area.]{IRI/LDEO:ClimateDataLibrary:2016}. In 2015, 167'377 OCV doses were distributed in the county of Juba during 6 days of a mass vaccination campaign started on July 31, 2015\cite[-3.2\baselineskip]{Abubakar:FirstUseGlobal:2015,Azman:EffectivenessOneDose:2016,Parker:AdaptingGlobalShortage:2017}, and are taken into account in the model.
%FFN: South Sudan has a lot of one dose
\newpage
\section{Cholera models}
\label{sec:meth}
\subsection{Generalized cholera model}
The proposed model is inspired by the model presented in the previous chapter: we have the susceptible $S$, infected $I$, and recovered $R$ compartments for individuals, with an additional variable $B$ describing the concentration of the bacteria in the environment. Previous modelling exercises had considered rainfall intensity $J(t)$ either to i) multiplicatively increase water \textit{contamination} with bacteria shed by infected individuals\cite{Bertuzzo:ProbabilityExtinctionHaiti:2016,Pasetto:RealtimeProjectionsCholera:2017}, or ii) assumed that the rainfall multiplicatively increases the \textit{exposure} to contaminated water\cite{Eisenberg:ExaminingRainfallCholera:2013}. A generalized formulation of these cholera-forced models, wherein both formulations are nested, is proposed. It enables the systematic comparison of the effect of rainfall through these two different transmission pathways.
The model described in \textsc{Chapter~1} was designed to deal with weekly data. Here, given the daily temporal resolution at which incidence data was available for the 2015 Juba epidemic, a compartment of exposed individuals $E$ is added to the model structure. It describes the incubation period: the lag between the time of exposure and the onset of the symptoms. Moreover, to account for the vaccination campaigns that were deployed in Juba in August 2015, four compartments ($V^S$, $V^E$, $V^I$, and $V^R$) are added to describe the dynamics of vaccinated individuals and their removal from the pool of susceptibles.
The proposed generalized cholera model is described in fig.~\ref{diagram}, and formulated as:
\marginnote[7\baselineskip]{
\begin{tabular}{ll}
\toprule
symbol & signification \\
\hline
$\alpha$ & cholera induced mortality\\
$\mu$ & natural mortality\\
$\sigma$ &symptomatic fraction\\
$\gamma$ & infectious period \\
$\rho$ & duration of immunity \\
$\rho_v$ & duration of vaccine protection \\
$\eta$ & vaccine efficacy \\
$\theta$ & shedding rate \\
$\mu_B$ & bacteria decay in water \\
$J(t)$ & rainfall \\
\bottomrule
\end{tabular}\\Reminder of the parameters already described in \textsc{Chapter~1} (p. \pageref{eq:I2}).}
\begingroup
\allowdisplaybreaks
\begin{eqnarray} \label{eq:fullmodel}
S(t) &=& H - I(t) - E(t) - R(t) - V^S(t) - V^E(t) - V^I(t) - V^R(t) \label{eq:S2j} \\
\frac{dE}{dt} &=& \sigma F(t) S - (\phi + \mu +\nu) E \label{eq:E2j}\\
\frac{dI}{dt} &=& \phi E - (\gamma + \mu + \alpha) I \label{eq:I2j}\\
\frac{dR}{dt} &=& (1-\sigma) F(t) S + \gamma I - (\rho + \mu+\nu) R \label{eq:R2j}\\
\frac{dB}{dt} &=& - \mu_B B +\theta\left[1 + f_{\mathcal{c}}\left(J(t)\right) \right] (I+V^I)\label{eq:B2}\\
\frac{dV^S}{dt} &=& \nu S - \mu V^S+ \rho_{v} V^R - (1-\eta) F(t) V^S \label{eq:VS2j}\\
\frac{dV^E}{dt} &=& \nu E + \sigma (1-\eta) F(t) V^S-(\phi + \mu) V^E \label{eq:VE2}\\
\frac{dV^I}{dt} &=& \phi V^E -(\gamma + \alpha + \mu) V^I \label{eq:VI2j}\\
\frac{dV^R}{dt} &=& \nu R -(\mu +\rho_{v})V^R +\gamma V^I +(1-\sigma) (1-\eta) F(t) V^S\label{eq:VR2j}\; ,
\end{eqnarray}
\endgroup
where $F(t)$ takes into account both human-to-human transmission and nonlinear water-to-human transmission:
\begin{equation}
F(t) = \beta_B \underbrace{ \left[\frac{B}{K + B} \right] \bigg(1+f_{\mathcal{e}}\left(J(t)\right)\bigg)}_{\text{indirect transmission}} + ~\beta_{I} \underbrace{\frac{(I+V^I)}{H}}_{\text{direct transmission}}.
\label{eq:force2}
\end{equation}
The notation is consistent with \textsc{Chapter~1} and we refer the reader to the description of the processes provided on p. \pageref{eq:I2}, and reminded in the margin. In addition to the previously described dynamics, exposed individuals become symptomatic infected at a rate $\phi$, which corresponds to an average incubation period of $1/\phi\approx$ 1.5 days\cite{Azman:IncubationPeriodCholera:2013}. As few reliable data on changes in Juba's population are available for the years of interest, the total population $H$ is assumed to be constant, as in \eqref{eq:S2j}.
Moreover, the effect of rainfall is more precisely described. Functions $f_{\mathcal{c}}\left(J(t)\right)$ and $f_{\mathcal{e}}\left(J(t)\right)$ account for the rainfall effect respectively by increasing the bacteria contamination in the water reservoir or directly through amplifying the exposure in the force of infection.
\begin{figure}
\centering
\includegraphics{fig_cholera-rainfall/Lemaitre_ACTROP_2018_42_R1_fig1.png}
\caption[Transition diagram for the competing cholera models]{Transition diagram for the different cholera models considered, with the different variations \textsc{me}, \textsc{mr}, and \textsc{mec} indicated.}
\label{diagram}
\end{figure}
With the objective of assessing the importance of rainfall on cholera transmission, a generalization of the linear relation found in the litterature\footnote{Specifically the reference exposure \parencite{Eisenberg:ExaminingRainfallCholera:2013} and contamination \parencite{Rinaldo:Reassessment20102011:2012} models.} is proposed, by using a nonlinear function form for $f_{\mathcal{c,e}}\left(J(t)\right)$:
\begin{equation}
f_{\mathcal{c,e}}\left(J(t)\right)=\lambda_{\mathcal{c,e}} \left(\frac{J(t)}{\max_t J(t)}\right)^{\alpha_{\mathcal{c,e}}},
\label{eq:nonlinear_rain}
\end{equation}
\begin{marginfigure}
\centering
\includegraphics{fig_cholera-rainfall/alpha.png}
\margincaption[Non-linear rainfall effect on transmission]{$\left(\frac{J_t}{J_{max}}\right)^\alpha$ for different values of $\alpha$. For $\alpha = 1$, the linear model described in previous works is obtained, and larger $\alpha$ corresponds to an increased relative importance of heavy rainfall events.}
\label{fig:alpha}
\end{marginfigure}
where the subscripts $_\mathcal{c,e}$ respectively denote the effect of rainfall on exposure and contamination, $\max_t J(t)$ is the maximum recorded rainfall intensity during the epidemic, and the $\alpha_{\mathcal{c,e}}\ge0$ controls for the relative importance of different rainfall intensities in their effect on the force of infection. Indeed, since the ratio $\frac{J(t)}{\max_t J(t)} \in [0,1]$, for $\alpha_{\mathcal{c,e}} \gg 1$ the ratio will tend to $0$ for all small precipitation events, leaving only the effect of the strongest events, whereas for $\alpha_{\mathcal{c,e}} < 1$ all precipitation events will be assigned a similar weight in the force of infection (fig. \ref{fig:alpha}). The formulations found in litterature are recovered by setting $\alpha_{\mathcal{c,e}} = 1$. The flexibility allowed by eqn.~\eqref{eq:nonlinear_rain} allows to discriminate between rainfall effects along a continuum from acting on disease transmission regardless of intensity to a threshold-like effect for the largest events which could be associated to severe flooding causing damages to the city's water and sanitary system, for instance leading to sewer overflow.
\subsection{Competing transmission models}
The relevance of the two rainfall-driven transmission pathways is assessed here by comparing the following models:
\begin{itemize}
\item \textsc{mn} SIRB model without rainfall: $\lambda_{\mathcal{c}} = \lambda_{\mathcal{e}} = 0$, $\beta_{I} = 0$, as the null hypothesis for the importance of rainfall.
\item \textsc{mc} SIRB model where for rainfall enhances the \textit{contamination} of the water reservoir: $\lambda_{\mathcal{e}} = 0$, $\beta_{I} = 0$.
\item \textsc{me} SIRB model where rainfall increases the \textit{exposure} to bacteria: $\lambda_{\mathcal{c}} = 0$, $\beta_{I} = 0$.
\item \textsc{mec} SIRB model combining both approaches \textsc{mc} and \textsc{me} ($\beta_{I} = 0$). Both ways of accounting rainfall play a role simultaneously.
\end{itemize}
\marginnote[-8\baselineskip]{\textsc{me} models uses a similar transmission pathway as in \textcite{Eisenberg:ExaminingRainfallCholera:2013}, whereas \textsc{mc} models were described in \textcite{Rinaldo:Reassessment20102011:2012}. These models were adapted to the generalized configuration, \eg in eqn.~\eqref{eq:force2} precipitation enters in the term $\left(1+f_E \left(J(t)\right)\right)$ which entails water-to-human transmission also when $J(t)=0$, unlike the published version. The steps take to adapt these models are left in the postprint \textsc{si}, section 1.}
\begin{table}[h!]
\centering
\begin{tabular}{lccccc}
\toprule
Model & $\lambda_{\mathcal{c}}$ & $\alpha_\mathcal{c}$ & $\lambda_{\mathcal{e}}$ & $\alpha_\mathcal{e}$ & $\beta_I$ \\
\hline
\textsc{mn} & -- & -- & -- & -- & -- \\
\textsc{mnh} & -- & -- & -- & -- & $\times$ \\
\textsc{mc} & $\times$ & $\times$ & -- & -- & -- \\
\textsc{me} & -- & -- & $\times$ & $\times$ & -- \\
\textsc{mch}& $\times$ & $\times$ & -- & -- &$\times$ \\
\textsc{meh}& -- & -- & $\times$ & $\times$ &$\times$ \\
\textsc{mec}& $\times$ & $\times$ & $\times$ & $\times$ & --\\
\textsc{mech}& $\times$ & $\times$ & $\times$ & $\times$ & $\times$ \\
\bottomrule
\end{tabular}
\caption[Parameters considered in the eight compared models]{Parameters considered in the eight compared models. $\lambda$ and $\alpha$ characterize the functional forms considering the precipitation eqn.~\eqref{eq:nonlinear_rain}. $\beta_I$ is the exposure for human-to-human transmission. A dash -- indicate that the parameter is set to zero, whereas a cross $\times$ indicate that the parameter is used.}
\label{tab:models}
\end{table}
For each model is explored the possibility of adding explicitly direct, human-to-human transmission ($\beta_{I} > 0$), which is indicated with an \textsc{H} at the end of the model name: \textsc{mnh}, \textsc{mch}, \textsc{meh}, and \textsc{mech}. The different parameters associated with the considered models are summarized in tab.~\ref{tab:models}.
The eight models are compared on the basis of their ability to match the time series of daily reported cases during the cholera epidemic in Juba of 2015. The postprint \textsc{si}, tab. 1 summarizes which parameters are calibrated for each model and their prior distribution. The degrees of freedom of the models, $n_p$, vary from $n_p=7$ for \textsc{mn} to $n_p=12$ for \textsc{mech}. Given the low number of daily reported cases and their ensuing variability, a stochastic equivalent is also implemented\footnote{The sole difference in formulation is the introduction of overdispersion in the infection process. The force of infection is multiplied by a time-continuous white noise process, as detailed in the eqns. \eqref{overdisp1}--\eqref{overdisp2}.} of the deterministic ODE system, eqns. \eqref{eq:fullmodel}-\eqref{eq:VR2j}, formulated as a continuous-time partially observed Markov process model, accounting for both demographic and disease transmission stochasticities\cite{Breto:TimeSeriesAnalysis:2009}. The generalized stochastic model is structured exactly like the deterministic one and a similar stochastic model is presented in the next chapter. Hence stochastic model equations are left in the postprint \textsc{si}.% left in the \textsc{Appendix}, eqn. \eqref{eq:stochsys} on p. \pageref{eq:stochsys}.
\subsection{Models calibration}
\marginnote{For details about fixed and calibrated parameters, along with prior and posterior distributions, the reader is referred to the supplementary information of \fullcite{Lemaitre:RainfallDriverEpidemic:2019}.}
\paragraph{Initial conditions} The history of cholera epidemics in South Sudan, and particularly in Juba, plays an important role in the determination of the size of susceptible and recovered compartments at the beginning of the 2015 epidemic. These two compartments are also largely impacted by the rate of immunity loss $\rho$ of the recovered individuals and by the probability of asymptomatic infection $1-\sigma$\footnote[][2\baselineskip]{in the literature, values range from $\sigma=0.5$, meaning one asymptomatic per each symptomatic infected, to less than $\sigma=0.01$, corresponding to more than 99 asymptomatic infected per each symptomatic infected \parencite{Fung:CholeraTransmissionDynamic:2014}.}. The initial conditions must therefore be estimated for each parameter set considered during calibration.
To take into account the uncertainty associated with the immunity landscape of Juba, the initial number of recovered individuals in April 2014 is calibrated for each model. Then, the detailed daily data of suspected cases during 2014 is used to estimate the associated number of recovered individuals at the start of the 2015 epidemic. Suspected cases in 2014 are scaled to the reporting fraction and undergo an exponential decay with an average time of immunity loss $1/\rho$\footnote{Something similar has been done in \parencite{Pasetto:RealtimeProjectionsCholera:2017}.}. Simulations are then initialized on June 5, 2015 considering two exposed individuals, two infected, and the associated steady-state bacteria concentration.
\paragraph{Calibration} Calibration of the deterministic model is performed using a Markov Chain Monte Carlo based algorithm\footnote{Differential Evolution Adaptive Metropolis, \textsc{dream} \parencite{Vrugt:MarkovChainMonte:2016}, was developed to explore high-dimensional parameter spaces. Given the parameters' prior distribution, \textsc{dream} searches and selects new samples in the posterior distribution by using multiple \textsc{mcmc} chains that run in parallel and that jointly contribute to the computation of the proposal parameter samples. \textsc{Mcmc} chains converge toward the posterior probability distribution based on the log-likelihood function of the data given the model output.} , which draws samples from the posterior distribution of the parameters.
Inference on the stochastic model is performed using a frequentist iterated filtering algorithm\footnote{MIF \parencite{Ionides:InferenceDynamicLatent:2015} is a frequentist-based approach for identifying the maximum likelihood estimation. It has proved successful even for a range of complex models of cholera dynamics \parencite{King:InapparentInfectionsCholera:2008,Baracchini:SeasonalityCholeraDynamics:2017}. The MIF2 algorithm, which employs iterated Bayes maps, is implemented in the R package \texttt{pomp} \parencite{King:StatisticalInferencePartially:2015}. %For this exercise, it is configured using 120 initial parameter vectors for each model, built using Sobol sequences over the parameter bounds.
}. Both model were fit against the daily reported cases accounting for over- or under- reporting, and assuming a Poisson likelyhood\shortcite{Camacho:CholeraEpidemicYemen:2018}. Each datapoint at reporting date $t_i$, $y_{t_i}$, is assumed to belong to a Poisson distribution centered on the modeled incidence $C_{t_i}$\footnote{$C_{t_i}$ is computed as in eqn. \eqref{eq:C}, but for the $E \rightarrow I$ transition. For the deterministic model: $ C_{t_i} = \int_{t_i}^{t_{i+1}} \phi \left(E(t) + V^E(t)\right) dt$. For the stochastic model the reporting process is described %in eqn. \eqref{eq:stochrep}.
as $ C_{t_i} = [N_{EI}(t_{i+1}) - N_{EI}(t_i)] + [N_{V^EV^I}(t_{i+1}) - N_{V^EV^I}(t_i)] $, where $N_{AB}$ denotes the stochastic counting process of transitions between classes $A$ and $B$ .
} and its associated parameter vector $\boldsymbol{\theta}$ altered by the reporting rate $\epsilon$, as:
\begin{equation}
y_{t_i} \sim \text{Poisson}\left(\epsilon \,C_{t_i}(\boldsymbol{\theta})\right), \; \epsilon > 0,
\label{eq:obs}
\end{equation}
Models were then compared using the Bayesian Information Criterion (BIC), Bayes factors, and the likelihood ratio test for the nested models.
\section{Results}
\subsection{Model selection of rainfall effects on cholera transmission}
The summary statistics of the deterministic and stochastic models considered in the study are given in tab.~\ref{tab:stats}. Overall, the stochastic models outperform their deterministic counterparts for all model structures by $\approx 40$ log-likelihood units. Both model types agree on the significance of rainfall in explaining the time series of daily reported cases, in particular through the increased exposure pathway, although the specific ordering of the models differs between model types. The BFs for the deterministic models suggest a strong support for model \textsc{mec}, followed by model \textsc{me} ($BF^{-1}_{ME,MEH} = 0.16$), with very little support for all other models ($BF^{-1}_{\boldsymbol{\cdot},MEH}< 10^{-2}$ for all other models than ME). For the stochastic model, the BFs estimated with the BIC suggest the strongest support for model \textsc{me}, with the basic SIRB model coming in second with 5 times less evidence ($BF^{-1}_{MN,ME} \approx 0.15$).
\begin{table}[h!]
\centering
\small
\begin{tabular}{lcccccccc}
\toprule
& \multicolumn{4}{c}{Deterministic} & \multicolumn{4}{c}{Stochastic} \\
\cmidrule(rl){2-5}\cmidrule(rl){6-9}
Model & $n$ & $\hat{\ell}$ & BIC & ${BF}^{-1}$ & $n$ & \makecell{$\hat{\ell}$ \vspace{-.1cm} \\ \scriptsize{(s.e.)}} & BIC & $BF^{-1}$ \\
\cmidrule(rl){2-5}\cmidrule(rl){6-9}
\textsc{mn} & $7$ & $-368.62$ & $770.27$ & $3.1\times 10^{-5}$ & $8$ & \makecell{$-326.45$ \vspace{-.1cm} \\ \scriptsize{$(0.105)$}} & $690.65$ & $1.5\times 10^{-1}$ \\
\textsc{mnh} & $8$ & $-368.95$ & $775.64$ & $1.1\times 10^{-9}$ & $9$ & \makecell{$-327.52$\vspace{-.1cm} \\ \scriptsize{$(0.052)$}} & $697.51$ & $4.7\times 10^{-3}$ \\
\textsc{mc} & $9$ & $-358.32$ & $759.11$ & $5.5\times 10^{-3}$ & $10$ & \makecell{$-323.50$\vspace{-.1cm} \\ \scriptsize{$(0.037)$}} & $696.01$ & $2.5\times 10^{-2}$ \\
\textsc{mch} & $10$ & $-359.06$ & $765.30$ & $1.7\times 10^{-4}$ & $11$ & \makecell{$-324.89$\vspace{-.1cm} \\ \scriptsize{$(0.041)$}} & $701.68$ & $5.9\times 10^{-4}$ \\
\textsc{me} & $9$ & $-356.96$ & $756.40$ & $1.6\times 10^{-1}$ & $10$ & \makecell{$-319.81$\vspace{-.1cm} \\ \scriptsize{$(0.035)$}} & $687.38$ & $1$ \\
\textsc{meh} & $10$ & $-358.06$ & $763.30$ & $6.3\times 10^{-4}$ & $11$ & \makecell{$-320.64$\vspace{-.1cm} \\ \scriptsize{$(0.030)$}} & $693.18$ & $4.1\times 10^{-2}$ \\
\textsc{mec} & $11$ & $-356.87$ & $765.64$ & $1$ & $12$ & \makecell{$-320.17$\vspace{-.1cm} \\ \scriptsize{$(0.031)$}} & $696.96$ & $6.2\times 10^{-3}$ \\
\textsc{mech} & $12$ & $-357.55$ & $771.73$ & $2.4\times 10^{-6}$ & $13$ & \makecell{$-320.38$\vspace{-.1cm} \\ \scriptsize{$(0.024)$}} & $702.09$ & $4.8\times 10^{-4}$ \\
\bottomrule
\end{tabular}
\caption[Model comparison statistics]{Model comparison statistics. For each model is reported its number of parameters $n$, the associated estimated log-likelihood $\hat{\ell}$ (and its Monte Carlo standard error for the stochastic model), and the inverse of the Bayes Factor ($BF^{-1}$) with respect to the model with the largest evidence. The BFs for the deterministic models were computed directly from the parameters posteriors, whereas for the stochastic models they were estimated with the Bayesian Information Criterion (BIC) as $BF_{i} \approx e^{\frac{1}{2} \left( BIC_i - BIC_{min}\right)}$. The BIC for the deterministic models was computed using the maximum log-likelihood value visited with the MCMC algorithm across chains.}\label{tab:stats}
\end{table}
When considering only the BIC, model \textsc{me} ranks first for both the deterministic and the stochastic formulations. Interestingly, all models that include human-to-human transmission present smaller or equal log-likelihoods than their counterparts with only the bacteria compartment, which suggests that the data does not support both environmental and human-to-human transmission within the set of the models considered here.
The results of the nested LR-tests confirm the statistical significance of including rainfall in the cholera transmission models, with the effect on exposure better supported by the data in both model types than the effect on contamination. In the deterministic case, the extension of the basic SIRB (model \textsc{mn}) with rainfall effects were significant for all direct comparisons (fig.~\ref{fig:lltests},\textsc{a}). The addition of human-to-human transmission was not significant mostly due to the above-mentioned lower estimate of the log-likelihood in these models. When considering only a single effect of rainfall (either increasing exposure or contamination), \textsc{me} outperforms \textsc{mc} in terms of likelihood for the same number of parameters. Interestingly, the inclusion of rainfall-induced contamination in model \textsc{me} is rejected due to the very limited increase of the estimated log-likelihood of \textsc{mec}, in contrast with the BFs favoring the latter. Model \textsc{me} is thus the one retained by the LR-tests in the deterministic set of models. In the case of the stochastic models, the LR-tests also highlight the importance of the effect of rainfall on exposure rather than on contamination (fig.~\ref{fig:lltests},\textsc{b}). In fact, the much stronger performance of \textsc{mn} in comparison with its deterministic counterpart relative to all other model structures imposes a stronger condition for retaining additional transmission processes. Indeed, both models \textsc{mc} and \textsc{mch} were rejected when compared to \textsc{mn}, thus only models with rainfall-driven exposure were retained. As in the deterministic case, model \textsc{me} is the one finally retained due to the lack of significance of the inclusion of additional transmission processes. Note that the conclusion based on the LR-tests for the deterministic models should be taken with caution because the \textsc{mcmc} algorithm used for calibration does not directly aim at maximizing the likelihood, but rather at sampling from the posterior distribution of the parameters given the data and the model. Moreover, the best likelihood visited by the chains when sampling the posteriors that are used in the LR-tests is not a formal estimate of the models' likelihood. However, the fact that the LR-tests applied to both model types agree with the selection of \textsc{me} supports their use in both cases.
\begin{figure*}
\centering
\includegraphics[width = \textwidth, trim = 11mm 19mm 10mm 12mm, clip]{fig_cholera-rainfall/Lemaitre_ACTROP_2018_42_R1_fig3.png}
\caption[Likelihood ratio tests of model nesting]{Likelihood ratio tests of model nesting. The LL-tests were computed for each nested pair of models $\{\mathcal{M}_0, \mathcal{M}_1\}$, with parameter vectors $\boldsymbol{\theta}^0,\boldsymbol{\theta}^1$, for which at least on of the parameters that is not null is $\boldsymbol{\theta}^1$ is equal to $0$ in $\boldsymbol{\theta}^0$. Each model is labeled with its associated estimated maximum log-likelihood value, $\hat{\ell}$, for the deterministic (A) and the stochastic (B) models, and linked based on whether the likelihood ratio is significantly (full black lines) or not (dashed gray lines) at the 5\% level. The absence of lines indicates a lower $\hat{\ell}$ for the more complex model.}
\label{fig:lltests}
\end{figure*}
Both statistical methods for model comparison therefore agree about the importance of the effect of intra-seasonal rainfall on the exposure to transmission during the 2015 cholera epidemic in Juba. For the deterministic type of models, the BFs suggest stronger support for model \textsc{mec}, and the LR-tests for \textsc{me}, whereas for the stochastic models both the BIC-based estimates of BFs and the LR-tests favor \textsc{me}.
\subsection{Intra-seasonal rainfall events and the 2015 Juba epidemic}
The comparison between the estimated output cases computed by the basic SIRB model (\textsc{mn}) and the most significant rainfall-based processes (\textsc{mec} and \textsc{me} for the deterministic and stochastic types, respectively) highlight the importance of rainfall in retrieving the second epidemiological peak (fig.~\ref{fig:sim}). Both deterministic and stochastic models fit well the general trend of the data, but they underestimate the large number of reported cases on July 19 (65 cases). Instead, the more complex models \textsc{me} and \textsc{mec} follow the SIRB dynamics and then are forced by the precipitation that occurred on July 18 (33 mm/d) toward the epidemiological peak.
%
Model calibration results suggest that precipitations with smaller intensities did not have a strong impact on cholera transmission during the 2015 epidemic in Juba. The exponents $\alpha_{\mathcal{c}}$ and $\alpha_{\mathcal{e}}$ were found to be systematically larger \mbox{than 1.} % (as shown by posteriors of the deterministic models in fig. S.2 and the Monte Carlo confidence intervals for the stochastic \textsc{me} in fig. S.4 of the SI) --> appendix or not ? DECISION
Thus, in the considered epidemic, the nonlinear function used to account for rainfall in the model, eqn.~\eqref{eq:nonlinear_rain}, helps to isolate the contribution of the largest rainfall.
The best measures of fit computed for the stochastic \textsc{me} (see tab.~\ref{tab:stats}) are thus explained by a larger sensitivity to precipitation, which causes the match between the mean of the simulated cases and the data during the second peak.
\begin{figure}[ht]
\centering
\includegraphics{fig_cholera-rainfall/Lemaitre_ACTROP_2018_42_R1_fig4.png}
\caption[Fit of the different models]{Simulations of the \textsc{mn} (A-B), \textsc{me} (C-D) and \textsc{mec} (E-F) models. Simulations for the deterministic versions (A,C,E) are given by the mean (blue dashed line), median (blue full line) and 95\% simulation envelops (blue ribbon) of 100 simulations of the measurement model for each trajectory from 100 samples from the posteriors of model parameters against reported daily cholera cases (red line and dots). Simulations from the stochastic models (B, D, F) are given for 10'000 simulations of the stochastic process and measurement models.}
\label{fig:sim}
\end{figure}
Comparing the two model types, stochastic results have a larger 95\% confidence interval, which better encompasses most of the data. In particular, both epidemiological peaks are well captured by the stochastic models, while the deterministic results systematically underestimate them. Two factors contribute to this result: the intrinsic stochastic nature of the model, which requires the simulation of various model runs for the same set of parameters, and the noise that necessarily perturbs the force of the infection yielding an overdispersion in infections. The standard deviation of such (assumed white) noise is estimated in each stochastic model, and it is interesting to note that the MLE obtained for \textsc{ME} is slightly smaller than in \textsc{MN} (0.028 versus 0.022), highlighting again that the data are retrieved with a lower uncertainty when rainfall is included in the model. This is evident in fig.~\ref{fig:sim}, where the width of the 95\% confidence interval of models \textsc{me} and \textsc{mec} is smaller with respect to \textsc{mn}.
Finally, despite having different BFs, the deterministic models \textsc{me} and \textsc{mec} are qualitatively similar in terms of output response, indicating that the recorded changes in the log-likelihood function do not correspond to qualitative changes in the output.
\section{Discussion}
\label{sec:disc}
In this chapter, a general mechanistic SIRB-based epidemiological model is developed to evaluate the relevance of rainfall in the amplification of cholera transmission, focusing on the 2015 Juba outbreak. Two rainfall-based transmission processes were compared: the direct increase of the exposure to the contaminated water (model~\textsc{me})\cite{Eisenberg:ExaminingRainfallCholera:2013}, and the increase of water contamination by flooded open defecation sites (model~\textsc{mc})\cite{Rinaldo:Reassessment20102011:2012}. In addition, human-to-human transmission is considered (models' name with \textsc{h}).
Regarding the epidemiological model, this study introduced two innovations with respect to previous modeling attempts of cholera epidemics\cite{Bertuzzo:ProbabilityExtinctionHaiti:2016,Pasetto:RealtimeProjectionsCholera:2017}. First, the focus on daily incidence data, as opposed to weekly epidemiological reports commonly used in cholera studies, motivated the introduction of a compartment of exposed individuals, eqn.~\eqref{eq:E2j}, to account for the incubation period of the disease and, thus, the lag between the possibly rainfall-driven infection process and the manifestation of the symptoms resulting in the time-series of daily reported cases at our disposal. %Such compartment had not been considered in previous cholera modelling works because they were most often working with weekly case reports. It was deemed necessary in this study to properly align the daily rainfall forcing with the daily reported cases.
Second, a non-linear version of rainfall driver, in the form of a power-law controlled by a single parameter, was introduced to generalize the previous linear dependence. Such formulation has the flexibility to either emphasize the impact of the largest rainfall events or give equal weight to all non-zero rainfall intensities. %Suffice here to mention that the added parameter is properly discounted in the formal model comparison carried out here.
%The idea is that large rainfall events may have a larger role in the deterioration of the sanitation conditions and, thus, in increasing the risk of transmission. The power low used for this goal allows, through its parameter $\alpha$, to emphasize large precipitation events or to flatten the effect for all precipitation.
All model assumptions were compared for both deterministic and stochastic model types, to draw more general conclusions. The statistics and tests used to compare the model results (tab.~\ref{tab:stats} and fig.~\ref{fig:lltests}) supported the significance of rainfall effects during the 2015 epidemic in Juba. In fact, results showed that for both model types there exists a significant positive effect of including rainfall drivers, in particular because standard SIRB models were not able to reproduce the second epidemiological peak of reported cases that occurred in July during the recession period. All models considering rainfall, instead, showed an increase of the number of cases in correspondence of the second peak, which was due to the large rainfall rates that occurred during the previous days (fig.~\ref{fig:sim}).
This difference in the simulated responses of models considering (or not) rainfall lead to stronger support for rainfall-based models. Due to the small variations among the likelihoods of rainfall-based models, however, it is not straightforward to conclude the best way to include the rainfall effect (tab.~\ref{tab:stats}).
Models with the minimum BIC were those considering the increase in exposure (model \textsc{me}) for both the stochastic and the deterministic model types. For the deterministic models, the computation of the Bayes Factors (BFs), which should provide a direct estimation of the model probability, suggested the selection of the model combining exposure and contamination processes (model \textsc{mec}). However, this information criterion might be unstable due to numerical issues and oscillations in the \textsc{mcmc} used for calibration\cite{Raftery:EstimatingIntegratedLikelihood:2007}. Because the models' outputs were similar for \textsc{mec} and \textsc{me} (fig. \ref{fig:sim}), it is advised to select the approach endowed with less parameters\footnote{in an attempt to maximize the model predictive accuracy we want to reduce potential overfitting.}, in this case~\textsc{me}, as indicated by the BIC. Note that the inclusion of human-to-human transmission was not statistically relevant in this modeling exercise.
\marginnote{It is also possible that the stochastic model has benefitted from improved identifiability. It has been shown that stochastic models can often extract more information about parameters than deterministic counterparts, see:
\fullcite{Browning:IdentifiabilityAnalysisStochastic:2020}.}
The comparison between the likelihoods of the two models' types (deterministic and stochastic) showed that considering the stochasticity of the processes improves the model results (tab.~\ref{tab:stats}). This suggests that also deterministic models should include a stochastic term in the computation of the force of infection, eqn.~\eqref{eq:force2}, which might increase the flexibility of the outputs.
Several limitations should be considered when analyzing the present results. The calibration exercise attempted in this study considered daily rainfall and cholera reported cases, which are characterized by significant random fluctuations that might partly cloud the description of the underlying infection processes.
Small random delays in reporting could change the infection curve and thus the effect of rainfall. This issue was partially addressed by considering the exposed compartment, eqn. \eqref{eq:E2j}, for simulation of the incubation period, and unknown reporting rate $\epsilon$ for the observed cases.
Here, it is aimed to reproduce the epidemic by modeling epidemiological transmission processes. %As usual for modeling studies, the level of description of such processes is chosen by balancing model complexity and accuracy: an increase in the number of modeled processes results in a larger number of unknown parameters, and thus a more complex calibration of the model.
While non-linear rainfall effects and possible over-reporting are taken into account, human mobility effect\cite[-4\baselineskip]{Gatto:GeneralizedReproductionNumbers:2012,Bertuzzo:SpatiallyExplicitModels:2010,Mari:PredictiveAbilityMechanistic:2015,Perez-Saez:ClimatedrivenEndemicCholera:2017} is not considered, which could help modeling the importation of infected individual. Moreover, in this model asymptomatic infected individuals do not contribute transmission. From a modeling viewpoint, these unaccounted processes were compensated by the calibration procedure, at the loss of predictive power.
The prior bounds to be assigned to parameters are typically wide\cite{Akman:ExaminationModelsCholera:2016} because the rates governing transmission processes are highly dependent on the specific epidemiological context so that somewhat contradictory values had been estimated in the literature. These considerations, together with the intrinsic noise affecting recorded cases, underlie the possibility that some of the model parameters might be unidentifiable\cite{Eisenberg:IdentifiabilityEstimationMultiple:2013}, in the sense that different parameter combinations would yield the same model output (also called equifinality).
%The exploration of the posterior parameter distribution using an \textsc{mcmc} approach allowed us to evince the possible correlation among parameters that were well identified by the data, with the main risk of the algorithm getting trapped in a local minimum of the fitness landscape (the distribution of parameters). The posterior probability distributions of the model parameters (see the postprint \textsc{si}, section S4) are associated with the model uncertainty and were here explored by the chains of the \textsc{mcmc} calibration.
The lack of available data prevented us to include the effects of the overall efforts towards WaSH improvement in this study, but vaccination is implemented in order to eliminate one possible covariate of the rainfall effect.
Despite these limitations, the proposed model comparison using both a deterministic and a stochastic model gave coherent results. The agreement of the two modeling types strengthened the results regarding the importance of rainfall patterns to significantly affect the development of cholera cases in time.
Overall, the findings of the study are consistent with the lessons learned in South Sudan with most of the transmission starting with the onset of the rainy season. In 2016 and 2017, cases in the dry season were observed and associated to the overexploitation of scarce water resources by nomadic herdsmen. This suggests that, as already observed, a general assessment of the relationship between precipitation and general waterborne or water-based disease infections is far from obvious and surely case-dependent. It has been argued, for example, that in the domain of water-based parasitic infections rainfall could not only boost disease transmission but also reduce it substantially\cite[-12\baselineskip]{McCreesh:ChallengesPredictingEffects:2013}, e.g. by increasing water flow (which in turn decreases habitat suitability in water). Rainfall patterns may also drastically affect human activities related to water contacts, thus potentially altering exposure and transmission risk\cite{Lai:SpatialDistributionSchistosomiasis:2015}. To that end, a hydrology-driven assessment cannot ignore certain characteristics, in particular the ephemeral or permanent nature of the waterways fostering contacts among pathogens and hosts\cite{Perez-Saez:HydrologyDensityFeedbacks:2016}. Also, temporal fluctuations of rainfall patterns may be particularly important in determining the seasonality of transmission\cite{Bertuzzo:HydroclimatologyDualpeakAnnual:2012,Bertuzzo:PredictionSpatialEvolution:2011,McCreesh:PredictingEffectsClimate:2015,Perez-Saez:HydrologyDensityFeedbacks:2016}.
%Subsuming the results obtained by this major computational exploration focused on the analysis of Juba's 2015 epidemics, it is concluded that rainfall patterns are fundamental drivers for epidemic cholera models, whether deterministic or stochastic, not only to capture seasonal trends but also to describe short-term fluctuations in the number of reported cases.