Conditions for transient epidemics of waterborne disease in spatially explicit systems
Abstract
Waterborne diseases are a diverse family of infections transmitted through ingestion of—or contact with—water infested with pathogens. Outbreaks of waterborne infections often show well-defined spatial signatures that are typically linked to local eco-epidemiological conditions, water-mediated pathogen transport and human mobility. In this work, we apply a spatially explicit network model describing the transmission cycle of waterborne pathogens to determine invasion conditions in metacommunities endowed with a realistic spatial structure. Specifically, we aim to define conditions under which pathogens can temporarily colonize a set of human communities, thus triggering a transient epidemic outbreak. To that end, we apply generalized reactivity analysis, a recently developed methodological framework for the study of transient dynamics in ecological systems subject to external perturbations. The study of pathogen invasion is complemented by the detection of the spatial signatures associated with the perturbations to a disease-free system that are expected to be amplified the most over different time scales. Understanding the drivers of waterborne disease dynamics over time scales that are relevant to epidemic and/or endemic transmission is a crucial, cross-disciplinary challenge, as large portions of the developing world still struggle to cope with the burden of these infections.
1. Introduction
Waterborne diseases are infections caused by ingestion of (or, more generally, contact with) water contaminated by pathogenic organisms, ranging from micro- (protozoa, bacteria, viruses, algae) to macro-parasites (helminths like flatworms and roundworms). They still represent a major threat to human health, especially in low-income countries. Cholera and typhoid fever are among the best-known examples of potentially lethal waterborne diseases. Also, diarrhea, commonly associated with waterborne pathogens, is responsible for the deaths of about 1.5 million people every year, thus representing one of the leading causes of death, especially among infants and children in the developing world [1]. Unsafe water supply, lack of sanitation and poor hygienic conditions, which either directly or indirectly affect exposure and transmission rates, are crucial factors in determining the burden of waterborne infections [2].
Developing appropriate mechanistic tools is fundamental to understand, predict and control waterborne disease transmission. Typically, the development of such tools lies at the interface between mathematical biology, environmental science and epidemiology. A prototypical model for microparasitic waterborne infections was proposed by Codeço [3], who developed a system of three ordinary differential equations (ODEs) where, in addition to the compartments of susceptible (S) and infected (I) that characterize traditional models for microparasite transmission, one equation accounts for the population dynamics of pathogens (B) in the water reservoir used by the human community under study. Spatial dynamics were not accounted for in early studies on waterborne disease transmission. However, the spatial component of disease spread is at least as important as its temporal dynamics in determining local changes in the abundances of the aforementioned epidemiological compartments, especially in large-scale applications. The spatial spread of waterborne diseases is driven by both water-mediated transport [4] and human mobility [5]. The former typically occurs within a drainage basin, along natural or man-made network systems (i.e. rivers or canals, respectively), while the latter provides an efficient means of across-catchment (and possibly long-distance) pathogen dispersal. Therefore, spatially explicit modelling is crucial to understand how different spatial transmission mechanisms interact with each other and influence the spread of waterborne infections. It may also enable scientists and decision makers to reproduce with unprecedented accuracy real-world epidemiological dynamics (e.g. [6–10]), provide scenarios for future epidemic development over different time scales (e.g. [11–16]), and discuss the effectiveness of controls in real (e.g. [17–20]) and realistic (e.g. [21,22]) case studies. Spatial modelling also radically changed the traditional approach for the definition of conditions for pathogen establishment in complex applications, classically based on the evaluation of the basic reproduction number, R_{0}, defined as the average number of secondary infections caused by one infectious individual introduced in a completely susceptible community: in fact, it has been noted that values of R_{0} locally greater than one are neither necessary nor sufficient for outbreaks to occur when spatial interactions are a factor, as processes like hydrologic transport and human mobility can radically alter local transmission dynamics [23–27].
While the study of the asymptotic properties of spatially explicit transmission models can help design effective control strategies intended to permanently reduce pathogen spread (or even break disease transmission), understanding transient phenomena of potential epidemiological interest can be important as well, namely to possibly prevent transitory epidemics triggered by external perturbations to a system in which endemic transmission is not possible—i.e. a system in which the so-called disease-free equilibrium (DFE) is asymptotically stable. Neubert & Caswell [28] proposed a simple measure of a system’s short-term instability to small perturbations. Specifically, they defined reactivity as the maximum instantaneous rate at which perturbations to a stable steady state can be amplified. Although reactivity has been studied in several ecological contexts (see e.g. [29] and references therein), epidemiological applications are still relatively scarce: Hosack et al. [30] performed a complete reactivity analysis of Ross’s [31] malaria model; Chitnis et al. [32] used reactivity analysis to derive epidemicity thresholds for simple models of Rift Valley fever transmission; Woodall et al. [33] applied reactivity analysis to a host–pathogen system with culling of the host population; Mari et al. [29] studied short-term instabilities connected to disease transmission in spatially implicit metapopulations and Mari et al. [34] studied the reactivity properties of simple (spatially implicit) models for waterborne and water-related diseases.
In this work, we use the recently developed method of generalized reactivity (hereafter, simply, g-reactivity), which is specifically tailored for ecological and epidemiological systems [29]. The definition of g-reactivity is based on the analysis of the maximum amplification rate of external perturbations to a steady state as measured through a linear transformation of the state variables (the system output) that can be chosen so as to be epidemiologically relevant. This is achieved, in particular, by including all the infection-related components of the state space (such as the abundances of exposed or infected people, or the concentration of pathogens in an environmental reservoir) in the output, while at the same time excluding the others (such as the abundances of susceptible or recovered people; see [34]). Here, we extend the g-reactivity framework to encompass the analysis of waterborne disease dynamics in realistic, spatially explicit settings. Although reactivity has already been used in the analysis of spatially extended metapopulation models (e.g. for the study of Turing instabilities; see [35]), to the best of our knowledge the present paper represents one of the first attempts to study the reactivity properties of a realistic spatial system (see [36] for an application to marine protected areas), and the first to apply g-reactivity analysis in a spatially explicit context. Specifically, we aim to produce a complete characterization of the g-reactivity properties of the DFE of a model describing the transmission of waterborne microparasites in a metacommunity embedded in a realistic landscape. The analysis of transient dynamics is contrasted with the study of the asymptotic stability of the DFE (following [23,24]), to allow for a full comparison of the conditions leading to short- versus long-term waterborne pathogen invasion in spatially explicit systems.
The manuscript is organized as follows. In the next section, we outline a network model aimed to describe the transmission cycle of waterborne microparasites. Although general in its formulation, the processes described in the model are actually inspired by cholera transmission dynamics. In the following two sections, we derive spatially explicit conditions under which the DFE of the network system is g-reactive (so that suitable perturbations to the DFE can be temporarily amplified in an epidemiologically relevant system output, before eventually fading out) or unstable (so that pathogens can invade the community and establish permanently therein). Next, we study the geographical signatures of the largest-growing perturbations to the DFE, as well as the geography of epidemic outbreak. Outbreak and establishment conditions are then numerically analysed in spatially explicit systems, each consisting of a set of human settlements distributed along a realistic river network. Finally, a short discussion of the main results of the work, with a focus on eco-epidemiological implications, closes the paper.
2. The model
The study of waterborne pathogen transmission is here tackled by means of a network model [10] that has already been used as starting point for both theoretical studies on waterborne disease dynamics [5,24,26] and the analysis of real-world cholera epidemics [11–16,19,23]. The model describes local epidemiological, demographic and ecological processes, pathogen transport along water systems and the effects of short-term human mobility on disease propagation. Network nodes represent n human communities of assigned population size, arranged in a given spatial setting, and connected by hydrologic pathways and human mobility.
Let S_{i}(t) and I_{i}(t) be the local abundances of susceptible and infected individuals in each node i of the network at time t, and let B_{i}(t) be the concentration of pathogens in the local water reservoirs which human communities have access to. Epidemiological dynamics and pathogen transport over the hydrologic and human mobility networks can be described by the following set of 3n ODEs:
As for the human host population, the dynamics of the susceptible compartment in each community (first equation of model (2.1)) is described as a balance between population demography and infections due to exposure to the pathogen. The host population, if uninfected, is assumed to be at demographic equilibrium N_{i}, with μ being the human mortality rate. The parameter β_{i} represents the site-specific rate of exposure to contaminated water, and B_{i}/(K + B_{i}) is the dose–response function describing the probability of becoming infected due to the exposure to a concentration B_{i} of pathogens (with K being the half-saturation constant [3]). Exposure to contaminated water for susceptible people of community i can occur either in their home community (with probability $1-{m}_{i}^{S}$, with ${m}_{i}^{S}$ being the overall probability of exposure outside the home site i, as determined by the mobility of susceptible individuals) or elsewhere (with probability ${m}_{i}^{S}{Q}_{ij}$, with Q_{ij} representing the probability that water contacts taking place outside the home site i occur in site j ≠ i, Q_{ii} = 0). Other routes of infections, such as fast human-to-human transmission, which have been proposed for waterborne diseases like cholera, are here neglected for simplicity, but could be dealt with within the same modelling framework (e.g. [27]).
The evolution of the infected compartment (second equation of model (2.1)) is a balance between newly infected individuals and losses due to recovery or natural/pathogen-induced mortality, with δ and γ being the rates of disease-induced mortality and recovery from infection, respectively. Note that the recovered compartment is not modelled explicitly, so that individuals who recover from the acute phase of disease are simply removed from the population, as they were conferred life-long immunity to reinfection. For cholera, as an example, recent estimates place the duration of immunity in the range 2.3–3.0 years [13]. As such, loss of acquired immunity is unlikely to influence transient, short-term epidemic dynamics, which are the main focus of the present work. Although simplistic, the choice of neglecting the dynamics of recovered individuals is thus deemed reasonable for the problem at hand.
As for the pathogen population, the dynamics of the local concentrations of pathogens in the aquatic environment (third equation of model (2.1)) is given by a balance between water contamination, pathogen mortality and hydrologic transport. Pathogens are released in water (e.g. excreted) by infected individuals (from either the local community, with probability $1-{m}_{i}^{I}$, with ${m}_{i}^{I}$ being the overall probability of contamination outside the home site i, as determined by the mobility of infected individuals; or elsewhere, with probability ${m}_{\hspace{0.17em}j}^{I}{Q}_{\hspace{0.17em}ji}$) at a site-specific rate p_{i} and immediately diluted in a well-mixed local water reservoir of size W_{i}. Free-living pathogens are assumed to die at rate ν_{i}. They can also move between any two neighbouring nodes of the hydrologic network (say from i to j) at rate l_{i} and with probability P_{ij}.
Some of the parameters of model (2.1)—namely those related to human demography (μ) and the physiological response to the disease (δ, γ, K)—are assumed to be constant over the spatial scales considered in this study. All the other parameters are allowed to be possibly site-dependent. A summary of the state variables and parameters of model (2.1) is given in table 1. We finally note that setting ${m}_{i}^{S}={m}_{i}^{I}=0$ and l_{i} = 0 for all i’s in model (2.1) produces a set of n disconnected models accounting only for local disease transmission processes. Unless stated otherwise, all the analyses and results presented in the next sections refer to the full model accounting also for spatial coupling mechanisms.
symbol | definition |
---|---|
state variables | |
S_{i} | abundance of susceptible human hosts at site i (i = 1 … n) |
I_{i} | abundance of infected human hosts at site i |
B_{i} | concentration of pathogens in the water reservoirs of site i |
parameters: local processes | |
N_{i} | human population size in absence of disease at site i |
μ | baseline human mortality rate (site-independent) |
β_{i} | rate of exposure to contaminated water at site i |
K | half-saturation constant of dose-response function (site-independent) |
δ | disease-induced mortality rate (site-independent) |
γ | recovery rate (site-independent) |
p_{i} [θ_{i}] | rate of water contamination at site i [θ_{i} = p_{i}/K, rescaled contamination rate] |
W_{i} | size of the water reservoir of site i |
ν_{i} | pathogen mortality rate at site i |
parameters: spatial coupling mechanisms | |
${m}_{i}^{S}$ | overall probability of exposure outside home site i |
${m}_{i}^{I}$ | overall probability of contamination outside home site i |
Q_{ij} | probability of water contact at j conditional to occurring outside home site i |
l_{i} | hydrologic transport rate of pathogens at site i |
P_{ij} | probability of hydrologic pathogen transport between sites i and j |
parameters: output transformation | |
c_{I}_{i} | weight assigned to infected hosts at site i |
c_{B}_{i} | weight assigned to bacterial concentration at site i |
3. Conditions for short-term pathogen outbreak
To study short-term pathogen outbreaks, we seek conditions under which small perturbations to the DFE, the state of model (2.1) in which S_{i} = N_{i}, I_{i} = 0 and B_{i} = 0 for all i’s, can initially grow. If the DFE is unstable, that is if the generalized reproduction number ${\mathcal{R}}_{0}$ is larger than one (note that ${\mathcal{R}}_{0}$ can be worked out from the analysis of the Jacobian J_{0}, i.e. the matrix that describes the dynamics of the system linearized in a neighbourhood of the DFE; for details, see appendix A in the electronic supplementary material), virtually any perturbation can grow and generate an epidemic, which is eventually followed by the establishment of endemic pathogen transmission. More subtle is the case of a stable DFE (${\mathcal{R}}_{0}<1$), in which perturbations can either decay monotonically or undergo transient growth (thus possibly generating an epidemic wave) before eventually fading out. Neubert & Caswell [28] proposed a simple, yet extremely effective method to characterize the transient dynamics associated with a stable equilibrium of a linear (or linearized) system of ODEs after a pulse perturbation. Specifically, they defined as reactive those stable steady states for which there exist small perturbations that can yield a transient amplification of the Euclidean norm of the state vector. This definition has recently been extended to a generalized, fully anisotropic reactivity framework (g-reactivity, [29]).
Following Mari et al. [29], the DFE of model (2.1) is g-reactive if there exist small perturbations that can lead to a transient growth of the Euclidean norm of a suitable system output (y) that is linked to the full state of the system (x = [s^{T}, i^{T}, b^{T}]^{T}, a 3n-dimensional vector whose components s = [S_{1}, … , S_{n}]^{T}, i = [I_{1}, … , I_{n}]^{T} and b = [B_{i}/K, … , B_{n}/K]^{T} correspond to susceptible humans, infected humans and bacterial concentrations (properly rescaled); the superscript T indicates matrix transposition) through a linear transformation (y = Cx, where C is a full-rank q × 3n real matrix, q ≤ 3n; hence y is a q-dimensional vector). In other words, the DFE is g-reactive if
Using definition (3.1) and the output transformation matrix C from (3.2), it is possible to show (appendix B, electronic supplementary material) that the DFE of model (2.1) is g-reactive if
It is also worth noticing that the g-reactivity properties of the DFE can actually be evaluated based on a matrix of reduced order n, namely ${\mathbf{F}}_{\mathbf{0}}={\mathbf{T}}_{\mathbf{0}}+{\mathbf{E}}_{\mathbf{0}}+{\mathbf{E}}_{\mathbf{0}}^{\mathbf{S}}+{\mathbf{E}}_{\mathbf{0}}^{\mathbf{I}}+{\mathbf{E}}_{\mathbf{0}}^{\mathbf{S}\mathbf{I}}$ (details in appendix B), where
4. Identification of critical perturbations over different time scales
The analysis of the dominant eigenvalue of matrix H_{0} (or F_{0}) enables us to determine whether small perturbations to a stable DFE can be temporarily amplified in a suitable (i.e. epidemiologically relevant) system output. By the Rayleigh principle [37], the eigenvector associated with the dominant eigenvalue of H_{0} corresponds to the structure of the single perturbation (lying in the row space of C) characterized by the largest rate of amplification in the limit t → 0, i.e. the optimal perturbation at time 0 [28,38]. According to the Perron–Frobenius theorem for non-negative matrices [37], the dominant eigenvector of H_{0} is characterized by strictly positive components. Because model (2.1) is endowed with a well-defined spatial structure, the dominant eigenvector of matrix H_{0} can be interpreted as the geographic signature of the perturbation leading to the fastest outbreak in the short term. On the other hand, it is possible to show (see again appendix B) that the dominant eigenvector of F_{0} refers to the bacterial components of the state space, although only close to the transient epidemicity boundary ${\mathcal{E}}_{0}=1$.
One could also wonder whether there exist ways to determine the geographic signature of the perturbations leading to the largest epidemic amplification (measured as α_{y}(t) = ||y(t)||/||y(0)||) over longer time scales. In this respect, generalized stability theory (see again [38]) predicts that, in a linear (or linearized) system, the perturbation leading to the largest growth of the Euclidean norm of the system state (say, α_{x}(t) = ||x(t)||/||x(0)||) for t → ∞ can be obtained from the eigenvectors of the Jacobian matrix of the system, specifically as the conjugate of the biorthogonal of the leading eigenvector of the Jacobian. Going back to our problem, we know (appendix A) that J_{0} is block-triangular, with the infection-related components of model (2.1), namely infected people and free-living pathogens, basically constituting an isolated sub-system that is not influenced by the susceptible components of the model. The dynamics of this subsystem close to the DFE are thus completely determined by ${\mathbf{J}}_{\mathbf{0}}^{\mathrm{\prime}}$, that is the submatrix of the Jacobian J_{0} of the full system restricted to the infection-related state variables. Recalling that, with the definition of matrix C given in (3.2), susceptibles are assumed not to influence the system output, it can immediately be verified that ${\mathbf{H}}_{\mathbf{0}}=H({\mathbf{J}}_{\mathbf{0}}^{\mathrm{\prime}})=({\mathbf{J}}_{\mathbf{0}}^{\mathrm{\prime}}+{{\mathbf{J}}_{\mathbf{0}}^{\mathrm{\prime}}}^{\mathrm{T}})/2$ if c_{I} and c_{B} are both equal to the identity matrix of size n, U_{n}. In this case, in fact, α_{y}(t) corresponds to α_{x}(t) evaluated for the infection-related subsystem (in which x = [i^{T}, b^{T}]^{T}). Therefore, we can conclude that the largest-growing perturbation in the system output for t → ∞ is given by the conjugate of the biorthogonal of the leading eigenvector of ${\mathbf{J}}_{\mathbf{0}}^{\mathrm{\prime}}$, albeit only if c_{I} = c_{B} = U_{n}.
We finally remark that the asymptotic results drawn from generalized stability theory can provide meaningful indications about the behaviour of a nonlinear system around an equilibrium point only if the associated Jacobian can describe the dynamics of the system over relatively long time scales. In the problem at hand, this may prove true for a stable DFE, in particular for parameter combinations for which the initial depletion of the susceptible compartment is not too fast. If this is the case, then the asymptotic optimal perturbation based on the eigenvectors of matrix ${\mathbf{J}}_{\mathbf{0}}^{\mathrm{\prime}}$ may represent a heuristic upper bound to the long-term amplification of generic small perturbations to a stable DFE.
5. Outbreak and establishment conditions in a river network
5.1. Application to realistic settings
To analyse outbreak/establishment conditions in a realistic setting, we assume that human communities constitute the nodes of a so-called optimal channel network (OCN), i.e. a mathematical structure characterized by scaling forms that closely conform to the observed geomorphological features of real river networks [39–42]. Because of this feature, OCNs have often been used as a template for the structure of the landscape in epidemiological applications (e.g. [4,5,7,24,26,43]). We further assume, without loss of generality, that the OCN landscape is embedded in a square of unitary side (this condition can be relaxed by imposing periodic boundary conditions). We consider OCNs endowed with three different spatial configurations, characterized by the position of the outlet (figure 1). OCNs have been generated following the algorithm described in Bertuzzo et al. [44,45], based on the procedure proposed by Rigon et al. [46]. Since the generation of OCNs is an intrinsically stochastic process, we consider several (16) replicas for each network geometry. The total number of network nodes (n ≈ 500) is preserved in each geometry and replica.
As for water-mediated pathogen movement, we apply conservative transport everywhere except for the outlet, from which pathogens are removed from the network. The specification of matrix P = [P_{ij}] thus goes as follows. Let p^{d} be the probability of downstream transport (and 1 − p^{d} the probability of upstream transport). If h is a headwater node and j its downstream neighbour, P_{hj} = 1 (reflecting boundary); for the inner nodes i of the network, P_{ij} = p^{d} if j is the downstream neighbour, or ${P}_{ij}=(1-{p}^{d})/{n}_{i}^{u}$ if j is one of the ${n}_{i}^{u}$ upstream neighbours of node i, with ${n}_{i}^{u}=2$ for most i’s; at the outlet, o, pathogens are discharged from the network with probability p^{d} (absorbing boundary), while ${P}_{oj}=(1-{p}^{d})/{n}_{o}^{u}$ if j is one of the ${n}_{o}^{u}$ upstream neighbours of node o. Note that absorbing conditions prevent pathogens from accumulating at the network outlet.
On top of hydrologic connectivity, a second mobility layer accounting for human movement has also to be specified in model (2.1). To that end, pairwise movement probabilities Q_{ij} are described through a gravity model [47,48] in which the attractiveness of node j for node i is assumed to be directly proportional to the population size of j and inversely proportional to the distance d_{ij} between the two nodes (through an exponential kernel with scale factor D), i.e. Q_{ij} ∝ N_{j} exp (−d_{ij}/D) (if i ≠ j, Q_{ii} = 0). Movement probabilities Q_{ij} constitute the entries of the human mobility matrix (Q = [Q_{ij}]). To stipulate that Q is row-stochastic (i.e. a matrix in which rows sum up to one), outgoing mobility fluxes are normalized by $\sum _{k\ne i}^{n}{N}_{k}\mathrm{exp}(-{d}_{ik}/D)$. Different mobility models can be easily accommodated in the formalism of system (2.1), provided that human mobility may be expressed in terms of movement probability (as quantified by ${m}_{i}^{S}$ and ${m}_{i}^{I}$) and trip distribution (e.g. in the form of an origin-destination matrix, as quantified by Q).
5.2. Numerical analysis of threshold conditions
Figure 2a shows the g-reactivity and stability properties of the DFE of model (2.1) as a function of the human exposure and contamination rates (β and θ = p/K, respectively) when all model parameters, including the weights given to infected prevalence and bacterial abundance in the different communities, are assumed to be homogeneously distributed in space (therefore, N = NU_{n}, W = WU_{n}, β = βU_{n}, θ = θU_{n}, ν = νU_{n}, m^{S} = m^{S}U_{n}, m^{I} = m^{I}U_{n}, l = lU_{n}, c_{I} = c_{I}U_{n} and c_{B} = c_{B}U_{n} are all scalar matrices). Results refer to the leftmost OCN replica shown in figure 1a, but the g-reactivity and stability thresholds evaluated with all the other OCN structures are virtually indistinguishable from those of figure 2a. This shows that, for an assigned specification of hydrologic transport and human mobility, the g-reactivity and stability properties of model (2.1) may be quite robust to variations of the underlying spatial domain, concerning e.g. the fine-scale structure of the OCN landscape or the position of the outlet node.
If both the exposure rate β and the contamination rate θ are small in magnitude, the DFE is stable (${\mathcal{R}}_{0}<1$)—which prevents long-term pathogen establishment—and non-g-reactive (${\mathcal{E}}_{0}<1$, evaluated for c_{I} = 1 and c_{B} = 1)—which rules out even transient epidemic outbreaks in the community. The parameter regions characterized by either relatively large values of β and small values of θ, or relatively large values of θ and small values of β correspond instead to a stable, g-reactive DFE (${\mathcal{R}}_{0}<1$, ${\mathcal{E}}_{0}>1$). In this case, transient epidemic waves can be triggered by suitable perturbations to the DFE, provided that either exposure or contamination is sufficiently large, but long-term pathogen establishment and endemic transmission are not possible. For these to happen, both β and θ need to be sufficiently large.
We also note from figure 2a that higher values of β and θ are needed in the network model (which includes pathogen transport and human mobility) than in a spatially implicit setting (i.e. in a model with l_{i} = 0 and ${m}_{i}^{S}={m}_{i}^{I}=0$ for all i’s) for the DFE to be g-reactive/unstable. This is specifically due to the presence of an absorbing boundary for hydrologic pathogen transport, which makes this spatial coupling mechanism non-conservative. As a matter of fact, if hydrologic transport of pathogens is negligible—or otherwise dominated by human mobility—a different result can be found, namely that epidemic/endemic transmission requires lower values of β and θ in a spatially explicit (rather than in a spatially implicit) setting (e.g. [24]).
The definition of matrix C given in (3.2) clearly influences the g-reactivity properties of the DFE. In fact, choosing different values for c_{I} or c_{B} can imply a change in the classification of the DFE from g-reactive to non-g-reactive, or vice versa (figure 2b). Therefore, g-reactivity classification is not absolute: in order to be meaningful, it requires a suitable (i.e. epidemiologically motivated) design of the output transformation.
The role played by the transport/mobility parameters in triggering disease epidemicity or endemicity is shown in figure 3. High values of the hydrologic transport rate l and the downstream transport probability p^{d} are associated with a stable, non-g-reactive DFE, while transient epidemics and endemic transmission can be found for lower values of l and/or p^{d} (panel a). On the other hand, human mobility promotes both short-term outbreaks and long-term pathogen establishment (panel b). In fact, for low levels of human mobility (small m^{S} and m^{I}) the pathogen cannot invade the system and the DFE is non-g-reactive. Finally, it is possible to study the interplay between water-mediated pathogen transport and human mobility (panel c). Again, high values of the hydrologic transport rate are associated with a stable, possibly non-g-reactive DFE. Conversely, high human mobility can lead to a g-reactive or unstable DFE. This general picture remains qualitatively unchanged for different values of the baseline exposure and contamination rates, although quantitative details can obviously vary. Interestingly, when looking at the parameters concerning spatial coupling mechanisms, we find that network topology may indeed influence the g-reactivity and stability properties of the DFE of model (2.1), especially in parameter regions characterized by relatively high hydrologic transport and/or human mobility (electronic supplementary material, figure S1).
5.3. Amplification of small perturbations to the DFE over different timescales
To verify that the analysis of the dominant eigenvector of matrix H_{0} does indeed make it possible to identify perturbations that are critical in terms of their initial rate of amplification, we solve numerically system (2.1) with different initial conditions and different OCN configurations. Specifically, we consider small perturbations in the form
One could argue that in real epidemics most perturbations of the DFE are likely to appear in the form of localized imports of either infected humans or bacteria. Figure 4b shows simulations of model (2.1) where the response of the system to different perturbations is measured by means of the output amplification α_{y}(t). Interestingly, the optimal perturbation at time 0 represents an empirical upper bound for the amplification of point-source perturbations of a stable (yet g-reactive) DFE—not only in the short term (as theoretically expected), but also over longer time scales. Actually, for the parameter setting considered here (characterized, in particular, by a relatively high value of the exposure rate, β, and a relatively low value of the contamination rate, θ), randomly localized imports of infected human hosts yield transient responses that wane monotonically over time, whereas randomly localized exogenous pathogen loads can lead to a temporary amplification of the system output. These transient responses are topped by the response associated with the optimal perturbation at time 0 even over relatively long time scales (almost four months in this example). However, model simulations show that the optimal perturbation at time 0 is far from being the most amplified perturbation over longer time scales, actually being taken over by both random, distributed perturbations (e.g. in the form of vector w_{rand} defined above) and the asymptotic optimal perturbation. The latter, in particular, is indeed found to be among the most amplified perturbations for long time scales, over which linearization could indeed be expected to fail. All these results are robust to changes in the underlying OCN configurations (compare figure 4 with figure S2, in electronic supplementary material).
5.4. Analysis of spatial patterns
The spatial signatures of the time-0 optimal perturbation are shown in the left panels of figure 5 (infected components only, but note that bacterial components are qualitatively similar to those shown in the figure) for different parameter combinations, all leading to a stable, g-reactive DFE. In the first two of the four considered parameter sets (panels a and b; the two parameter combinations are those indicated in figure 2a), the backbone of the optimal perturbation at time 0 is localized in a neighbourhood of the network outlet—especially so in the second case, which is characterized by smaller exposure and larger contamination rates. On the contrary, for parameter settings identifying slower hydrologic pathogen transport (c) or higher levels of human mobility (d), the backbone of the optimal perturbation at time 0 covers larger regions of the landscape, stretching relatively far off upstream the network outlet (which still remains an important focal point in the geography of this critical perturbation), and away from the main course of the stream network. Qualitatively consistent results are found with different OCN configurations (electronic supplementary material, figures S3 and S4).
Interestingly, the spatial structure of the asymptotic optimal perturbation (middle panels of figure 5) seems to be quite robust to changes in the parametrization of hydrologic transport and human mobility, and suggests that the contamination of headwaters or, more generally, of reaches that lie far away from the network outlet may yield relatively large and long-lasting (albeit transient) epidemic waves. The rationale behind this result is inherently spatial in nature, and can be summarized as follows. In a homogeneous spatial setting, large outbreaks are expected for perturbations that develop into transitory epidemics covering large areal extents in their spatio-temporal evolution; in the absence of sustained transmission and in the presence of directional pathogen movement, this happens to be true for initial conditions that are concentrated farthest from the outlet, as they can indeed generate transient epidemics in which the peak values of local infected prevalence and pathogen concentration travel downstream as hydrologic transport move pathogens towards the absorbing outlet, thus progressively spanning the whole network.
Finally, the spatial structure of the dominant eigenvector of matrix ${\mathbf{J}}_{\mathbf{0}}^{\mathrm{\prime}}$ (right panels of figure 5), which describes the geography of the epidemic outbreak as it fades out after the transient phenomena associated with the initial perturbation to the DFE have vanished (appendix A; see also [23,24]), turns out to be pretty robust to changes in model parameters as well, and indicates again the network outlet as a focal point of the spatial epidemic pattern. Interestingly, the perturbation that optimally excites the dominant eigenvector of the Jacobian (the one that dominates in the long run, shown in the middle panels) may bear little resemblance to the eigenvector itself. In fact, as noted by Farrell & Ioannau [38], in generic, non-normal systems (i.e. systems in which the eigenvectors of the associated Jacobian matrix do not form an orthonormal basis), an eigenvector and its biorthogonal may differ greatly.
6. Discussion and conclusion
In this work, we have analysed short- and long-term invasion dynamics for waterborne pathogens in realistic riverine landscapes. Specifically, using a well-established, spatially explicit transmission model, we have provided conditions for the DFE of the system (the state in which neither infected human hosts nor free-living pathogens are present in any of the communities included in the region under study) to be g-reactive [29]. Perturbations to a stable, g-reactive DFE may grow in the system output, defined as a linear transformation of the system state that can be chosen to be epidemiologically relevant [34], before eventually fading out. From an epidemiological perspective, a g-reactive DFE may be associated with transient epidemic waves, while an unstable DFE will lead to sustained epidemic dynamics followed by the establishment of endemic transmission.
Understanding the conditions under which short- and/or long-term pathogen transmission is possible in an assigned spatial setting bears important implications for predicting—and possibly controlling—disease spread. Our spatially explicit methods could in fact be used to determine conditions for the occurrence of real-world outbreaks, as well as to assess which communities will be hit with more strength during a waterborne disease epidemic [23]. Results of this kind could thus assist in the planning of interventions aimed to minimize outbreak risk by reducing human exposure and/or contamination, as well as in the management of emergencies in the aftermath of an epidemic outbreak. In this respect, g-reactivity analysis could serve as an ideal companion to asymptotic stability analysis in the definition of a multicriterial decision support framework (e.g. [49]) in which possible trade-offs between short- and long-term objectives can be fully disclosed and clearly examined.
The methods described in this paper can also help evaluate the spatial structure of the perturbations to the DFE that are expected to grow the most in the short run (and, under suitable conditions, in the mid-to-long run as well). Therefore, results from g-reactivity analysis of epidemiological models applied to real case studies could serve as heuristics to define upper bounds of the growth of small perturbations to the DFE, thus providing informed estimates about the maximum size of an outbreak before it goes off. Suggestions from g-reactivity analysis could also be used to optimize the structure of disease surveillance networks (e.g. [50]), specifically to make sure that critical perturbations can be readily identified—and effective interventions set up in a timely manner. For instance, surveillance proved fundamental in the containment of the cholera outbreak that struck Haiti in October 2010 and that has rapidly become one of the largest waterborne disease epidemics of the recent past [51,52].
It has to be remarked, however, that predictions drawn from g-reactivity analysis may prove relevant just shortly after the start of an epidemic in the presence of an unstable DFE. In this case, in fact, linearization is indeed expected to fail rapidly. Needless to say, this scenario is critical from an epidemiological standpoint, because virtually any small perturbation to an unstable DFE can lead to a large outbreak. For waterborne infections, factors contributing to the instability of the DFE are high exposure and contamination risk (say, because of inadequate water provisioning and sanitation infrastructures), the availability of suitable water reservoirs for the thriving of local pathogen populations, and human mobility. History retrospectively shows that all these conditions were unfortunately met in post-earthquake Haiti, where a single contamination event (involving a tributary of the Artibonite River, the largest Haitian fluvial system) from an external source was in fact sufficient to cause a devastating epidemic (e.g. [53–56]). G-reactivity and stability analyses can help understand the role played by different factors in the definition of transmission dynamics over a range of time scales also in the presence of an unstable DFE, and can be used to increase awareness about the risk of potentially catastrophic outbreaks—a necessary condition to ensure that preventable disasters like the Haiti cholera outbreak do not happen again in the future.
As far as surveillance is concerned, our analysis highlights the neighbourhood of the river network outlet as a critical region for disease spread: on the one hand, in fact, it may constitute the spatial backbone of the optimal perturbation at time 0; on the other, it represents a hotbed for epidemic spread also over longer time scales. These results, together with historical data showing that most outbreaks of cholera (a prototypical example among microparasitic waterborne diseases) originated in coastal regions ([57] and references therein), indicate that the river network outlet should represent a focal point for a surveillance network for waterborne infections and, possibly, for the deployment of healthcare resources. On the other hand, epidemic surveillance and control should not be geographically limited to the outlet region, as the upstream sections of the river network seem to also play an important role, e.g. in the definition of the asymptotic optimal perturbation. Epidemic surveillance and forecasting, design of ex-ante and ex-post interventions, and resource allocation for transmission control are all crucial objectives of public health policies. Such goals could thus greatly benefit from spatially explicit mathematical modelling.
An important aspect of pathogen transmission that has been neglected in this work is the seasonal variability of model parameters, such as the human exposure and contamination rates, the average pathogen lifespan, or the water reservoir volume [58]. Such temporal variability can result in marked fluctuations of the force of infection, with remarkable consequences for long-term pathogen invasion, as already shown by previous work on the stability properties of the DFE of waterborne disease models in spatially explicit and time-periodic settings [26]. On the other hand, reactivity analysis has never been performed in spatially explicit, time-varying environments. Indeed, the study of the interplay between seasonality and reactivity has started only recently, specifically in spatially implicit ecological systems (e.g. [59]). Because of the local properties of g-reactivity, though, it is safe to predict that seasonal fluctuations of the transmission parameters must play an important—yet possibly nontrivial to unravel—role in the short-term response of the DFE of an epidemiological model to small, time-confined perturbations. A future extension of this work may thus usefully concern the g-reactivity analysis of a seasonal model for waterborne disease dynamics in the framework of Floquet theory (see again [26]).
Although tailored here for waterborne microparasites, spatially explicit tools for g-reactivity and asymptotic stability can be applied to a suite of diverse epidemiological problems. In this respect, the closest application would be the study of g-reactivity for waterborne macroparasitic infections (such as schistosomiasis), which may possibly share similar hydroclimatological and socioeconomic drivers in spite of the different underlying transmission mechanisms [60], and for which spatially explicit modelling approaches are becoming increasingly available (e.g. [25,61–64]). Moreover, future studies may also focus on spatially explicit models for diseases with different infection pathways, including human-to-human, airborne, faecal-oral, sexual and vector-borne. Properly informed with the spatial mechanisms relevant to each transmission route, g-reactivity and stability analyses could help ecologists and epidemiologists better understand how pathogen transport and human mobility interact with local transmission mechanisms in shaping short- and long-term invasion dynamics, with important consequences for the effectiveness of infectious disease control.
Data accessibility
This work makes no use of original data. Sections 3 and 4 in the manuscript provide detailed instructions to reproduce the results presented in §5, where the model described in §2 is applied to theoretical riverine landscapes produced following procedures that are fully described elsewhere (references in §5). All analyses have been performed using numerical algorithms that are routinely implemented in standard scientific computing software. Specifically, MATLAB^{TM} version R2017a has been used for this work. A sample MATLAB^{TM} implementation of the instructions needed to evaluate the stability and g-reactivity properties of model (2.1) with the output transformation defined in (3.2) is provided in appendix C (electronic supplementary material).
Authors' contributions
L.M. conceived and coordinated the study, performed the numerical analyses and drafted the first version of the manuscript. R.C. contributed to drafting the manuscript. E.B. participated in numerical analyses. A.R. and M.G. contributed to conceiving and coordinating the study. All authors participated in the design of the study and the revision of the manuscript, and gave final approval for publication.
Competing interests
The authors have no competing interests.
Funding
L.M., R.C. and M.G. acknowledge support from Politecnico di Milano. L.M. and R.C. were also supported by Politecnico di Milano through the Polisocial Award programme, project MASTR-SLS. A.R. acknowledges funding from the ERC Advanced Grant RINEC 22761 ‘River networks as ecological corridors for species, populations, and pathogens of waterborne disease’ and from the Swiss National Science Foundation project ‘Optimal control of intervention strategies for waterborne disease epidemics’ (SNSF grant no. 200021_172578).
Acknowledgements
The authors wish to thank three anonymous reviewers for their useful comments on the original version of the manuscript.
References
- 1. World Health Organization. 2014 Preventing diarrhoea through better water, sanitation and hygiene. Technical Report, World Health Organization, Geneva, Switzerland. Google Scholar
- 2. World Health Organization. 2014 Global health observatory. Technical Report, World Health Organization, Geneva, Switzerland. Google Scholar
- 3.
Codeço CT . 2001 Endemic and epidemic dynamics of cholera: the role of the aquatic reservoir. BMC Infect. Dis. 1, 1. (doi:10.1186/1471-2334-1-1) Crossref, PubMed, Web of Science, Google Scholar - 4.
Bertuzzo E, Casagrandi R, Gatto M, Rodriguez-Iturbe I, Rinaldo A . 2010 On spatially explicit models of cholera epidemics. J. R. Soc. Interface 7, 321-333. (doi:10.1098/rsif.2009.0204) Link, Web of Science, Google Scholar - 5.
Mari L, Bertuzzo E, Righetto L, Casagrandi R, Gatto M, Rodriguez-Iturbe I, Rinaldo A . 2012 On the role of human mobility in the spread of cholera epidemics: towards an epidemiological movement ecology. Ecohydrology 5, 531-540. (doi:10.1002/eco.262) Crossref, Web of Science, Google Scholar - 6.
Bertuzzo E, Mari L, Righetto L, Gatto M, Casagrandi R, Rodriguez-Iturbe I, Rinaldo A . 2011 Prediction of the spatial evolution and effects of control measures for the unfolding Haiti cholera outbreak. Geophys. Res. Lett. 38, L06403. (doi:10.1029/2011GL046823) Crossref, Web of Science, Google Scholar - 7.
Bertuzzo E, Mari L, Righetto L, Gatto M, Casagrandi R, Rodriguez-Iturbe I, Rinaldo A . 2012 Hydroclimatology of dual-peak annual cholera incidence: insights from a spatially explicit model. Geophys. Res. Lett. 39, L05403. (doi:10.1029/2011GL050723) Crossref, Web of Science, Google Scholar - 8.
Finger F, Genolet T, Mari L, Constantin De Magny G, Manga NM, Rinaldo A, Bertuzzo E . 2016 Mobile phone data highlights the role of mass gatherings in the spreading of cholera outbreaks. Proc. Natl Acad. Sci. USA 113, 6421-6426. (doi:10.1073/pnas.1522305113) Crossref, PubMed, Web of Science, Google Scholar - 9.
Finger F, Knox AC, Bertuzzo E, Mari L, Bompangue D, Gatto M, Rodriguez-Iturbe I, Rinaldo A . 2014 Cholera in the Lake Kivu region (DRC): integrating remote sensing and spatially explicit epidemiological modeling. Water Resour. Res. 50, 5624-5637. (doi:10.1002/wrcr.v50.7) Crossref, Web of Science, Google Scholar - 10.
Mari L, Bertuzzo E, Righetto L, Casagrandi R, Gatto M, Rodriguez-Iturbe I, Rinaldo A . 2012 Modelling cholera epidemics: the role of waterways, human mobility and sanitation. J. R. Soc. Interface 9, 376-388. (doi:10.1098/rsif.2011.0304) Link, Web of Science, Google Scholar - 11.
Bertuzzo E, Finger F, Mari L, Gatto M, Rinaldo A . 2016 On the probability of extinction of the Haiti cholera epidemic. Stoch. Environ. Res. Risk Assess. 30, 2043-2055. (doi:10.1007/s00477-014-0906-3) Crossref, Web of Science, Google Scholar - 12.
Mari L, Bertuzzo E, Finger F, Casagrandi R, Gatto M, Rinaldo A . 2015 On the predictive ability of mechanistic models for the Haitian cholera epidemic. J. R. Soc. Interface 12, 20140840. (doi:10.1098/rsif.2014.0840) Link, Web of Science, Google Scholar - 13.
Pasetto D et al. 2018 Near real-time forecasting for cholera decision making in Haiti after Hurricane Matthew. PLoS Comput. Biol. 14, e1006127. (doi:10.1371/journal.pcbi.1006127) Crossref, PubMed, Web of Science, Google Scholar - 14.
Pasetto D, Finger F, Rinaldo A, Bertuzzo E . 2017 Real-time projections of cholera outbreaks through data assimilation and rainfall forecasting. Adv. Water Res. 108, 345-356. (doi:10.1016/j.advwatres.2016.10.004) Crossref, Web of Science, Google Scholar - 15.
Righetto L, Bertuzzo E, Mari L, Schild E, Casagrandi R, Gatto M, Rodriguez-Iturbe I, Rinaldo A . 2013 Rainfall mediations in the spreading of epidemic cholera. Adv. Water Res. 60, 34-46. (doi:10.1016/j.advwatres.2013.07.006) Crossref, Web of Science, Google Scholar - 16.
Rinaldo A et al. 2012 Reassessment of the 2010–2011 Haiti cholera outbreak and rainfall-driven multiseason projections. Proc. Natl Acad. Sci. USA 109, 6602-6607. (doi:10.1073/pnas.1203333109) Crossref, PubMed, Web of Science, Google Scholar - 17.
Chao DL, Halloran ME, Longini IM . 2011 Vaccination strategies for epidemic cholera in Haiti with implications for the developing world. Proc. Natl Acad. Sci. USA 108, 7081-7085. (doi:10.1073/pnas.1102149108) Crossref, PubMed, Web of Science, Google Scholar - 18.
Finger F et al. 2018 The potential impact of case-area targeted interventions in response to cholera outbreaks: a modeling study. PLoS Med. 15, e1002509. (doi:10.1371/journal.pmed.1002509) Crossref, PubMed, Web of Science, Google Scholar - 19.
Kühn J, Finger F, Bertuzzo E, Borgeaud S, Gatto M, Rinaldo A, Blokesch M . 2014 Glucose- but not rice-based oral rehydration therapy enhances the production of virulence determinants in the human pathogen Vibrio cholerae. PLoS Negl. Trop. Dis. 8, e3347. (doi:10.1371/journal.pntd.0003347) Crossref, PubMed, Web of Science, Google Scholar - 20.
Tuite AL, Tien J, Eisenberg M, Earns DJD, Ma J, Fisman DN . 2011 Cholera epidemic in Haiti, 2010: using a transmission model to explain spatial spread of disease and identify optimal control interventions. Ann. Intern. Med. 154, 593-601. (doi:10.7326/0003-4819-154-9-201105030-00334) Crossref, PubMed, Web of Science, Google Scholar - 21.
Gaythorpe K, Adams B . 2016 Disease and disaster: optimal deployment of epidemic control facilities in a spatially heterogeneous population with changing behaviour. J. Theor. Biol. 397, 169-178. (doi:10.1016/j.jtbi.2016.03.006) Crossref, PubMed, Web of Science, Google Scholar - 22.
Kelly MR, Tien JH, Eisenberg MC, Lenhart S . 2016 The impact of spatial arrangements on epidemic disease dynamics and intervention strategies. J. Biol. Dyn. 10, 222-249. (doi:10.1080/17513758.2016.1156172) Crossref, PubMed, Web of Science, Google Scholar - 23.
Gatto M, Mari L, Bertuzzo E, Casagrandi R, Righetto L, Rodriguez-Iturbe I, Rinaldo A . 2012 Generalized reproduction numbers and the prediction of patterns in waterborne disease. Proc. Natl Acad. Sci. USA 48, 19 703-19 708. (doi:10.1073/pnas.1217567109) Crossref, Web of Science, Google Scholar - 24.
Gatto M, Mari L, Bertuzzo E, Casagrandi R, Righetto L, Rodriguez-Iturbe I, Rinaldo A . 2013 Spatially explicit conditions for waterborne pathogen invasion. Am. Nat. 182, 328-346. (doi:10.1086/671258) Crossref, PubMed, Web of Science, Google Scholar - 25.
Gurarie D, Seto EYW . 2009 Connectivity sustains disease transmission in environments with low potential for endemicity: modelling schistosomiasis with hydrologic and social connectivities. J. R. Soc. Interface 6, 495-508. (doi:10.1098/rsif.2008.0265) Link, Web of Science, Google Scholar - 26.
Mari L, Casagrandi R, Bertuzzo E, Rinaldo A, Gatto M . 2014 Floquet theory for seasonal environmental forcing of spatially-explicit waterborne epidemics. Theor. Ecol. 7, 351-365. (doi:10.1007/s12080-014-0223-y) Crossref, Web of Science, Google Scholar - 27.
Tien JH, Shuai Z, Eisenberg MC, van den Driessche P . 2015 Disease invasion on community networks with environmental pathogen movement. J. Math. Biol. 70, 1065-1092. (doi:10.1007/s00285-014-0791-x) Crossref, PubMed, Web of Science, Google Scholar - 28.
Neubert MG, Caswell H . 1997 Alternatives to resilience for measuring the responses of ecological systems to perturbations. Ecology 78, 653-665. (doi:10.1890/0012-9658(1997)078[0653:ATRFMT]2.0.CO;2) Crossref, Web of Science, Google Scholar - 29.
Mari L, Casagrandi R, Rinaldo A, Gatto M . 2017 A generalized definition of reactivity for ecological systems and the problem of transient species dynamics. Methods Ecol. Evol. 8, 1574-1584. (doi:10.1111/2041-210X.12805) Crossref, Web of Science, Google Scholar - 30.
Hosack GR, Rossignol PA, van den Driessche P . 2008 The control of vector-borne disease epidemics. J. Theor. Biol. 255, 16-25. (doi:10.1016/j.jtbi.2008.07.033) Crossref, PubMed, Web of Science, Google Scholar - 31.
- 32.
Chitnis N, Hyman JM, Manore CA . 2013 Modelling vertical transmission in vector-borne diseases with applications to Rift Valley fever. J. Biol. Dyn. 7, 11-40. (doi:10.1080/17513758.2012.733427) Crossref, PubMed, Web of Science, Google Scholar - 33.
Woodall H, Bullock JM, White SM . 2014 Modelling the harvest of an insect pathogen. Ecol. Model 287, 16-26. (doi:10.1016/j.ecolmodel.2014.04.026) Crossref, Web of Science, Google Scholar - 34.
Mari L, Casagrandi R, Rinaldo A, Gatto M . 2018 Epidemicity thresholds for water-borne and water-related diseases. J. Theor. Biol. 447, 126-138. (doi:10.1016/j.jtbi.2018.03.024) Crossref, PubMed, Web of Science, Google Scholar - 35.
Neubert MG, Caswell H, Murray JD . 2002 Transient dynamics and pattern formation: reactivity is necessary for Turing instabilities. Math. Biosci. 175, 1-11. (doi:10.1016/S0025-5564(01)00087-6) Crossref, PubMed, Web of Science, Google Scholar - 36.
Aiken CM, Navarrete SA . 2011 Environmental fluctuations and asymmetrical dispersal: generalized stability theory for studying metapopulation persistence and marine protected areas. Mar. Ecol. Prog. Ser. 428, 77-88. (doi:10.3354/meps09079) Crossref, Web of Science, Google Scholar - 37.
Horn RA, Johnson CR . 1990 Matrix analysis. Cambridge, UK: Cambridge University Press. Google Scholar - 38.
Farrell BF, Ioannau PJ . 1996 Generalized stability theory. Part I: autonomous operators. J. Atmos. Sci. 53, 2025-2040. (doi:10.1175/1520-0469(1996)053<2025:GSTPIA>2.0.CO;2) Crossref, Web of Science, Google Scholar - 39.
Balister P, Balogh J, Bertuzzo E, Bollobás B, Caldarelli G, Maritan A, Mastrandrea R, Morris R, Rinaldo A . 2018 River landscapes and optimal channel networks. Proc. Natl Acad. Sci. USA, 115, 6548-6553. (doi:10.1073/pnas.1804484115) Crossref, PubMed, Web of Science, Google Scholar - 40.
Rinaldo A, Rigon R, Banavar JR, Maritan A, Rodriguez-Iturbe I . 2014 Evolution and selection of river networks: statics, dynamics, and complexity. Proc. Natl Acad. Sci. USA 111, 2417-2424. (doi:10.1073/pnas.1322700111) Crossref, PubMed, Web of Science, Google Scholar - 41.
Rinaldo A, Rodriguez-Iturbe I, Rigon R, Bras R, Ijjasz–Vasquez E, Marani A . 1992 Minimum energy and fractal structures of drainage networks. Water Resour. Res. 28, 2183-2195. (doi:10.1029/92WR00801) Crossref, Web of Science, Google Scholar - 42.
Rodriguez-Iturbe I, Rinaldo A, Rigon R, Bras R, Ijjasz–Vasquez E, Marani A . 1992 Fractal structures as least energy patterns – the case of river networks. Geophys. Res. Lett. 19, 889-892. (doi:10.1029/92GL00938) Crossref, Web of Science, Google Scholar - 43.
Carraro L, Mari L, Gatto M, Rinaldo A, Bertuzzo E . 2018 Spread of proliferative kidney disease in fish along stream networks: a spatial metacommunity framework. Freshw. Biol. 63, 114-127. (doi:10.1111/fwb.2018.63.issue-1) Crossref, Web of Science, Google Scholar - 44.
Bertuzzo E, Carrara F, Mari L, Altermatt F, Rodriguez-Iturbe I, Rinaldo A . 2016 Geomorphic controls on elevational gradients of species richness. Proc. Natl Acad. Sci. USA 113, 1737-1742. (doi:10.1073/pnas.1518922113) Crossref, PubMed, Web of Science, Google Scholar - 45.
Bertuzzo E, Rodriguez-Iturbe I, Rinaldo A . 2015 Metapopulation capacity of evolving fluvial landscapes. Water Resour. Res. 51, 2696-2706. (doi:10.1002/2015WR016946) Crossref, Web of Science, Google Scholar - 46.
Rigon R, Rinaldo A, Rodriguez-Iturbe I . 1994 On landscape self-organization. J. Geophys. Res. 99, 11 971-11 993. (doi:10.1029/93JB03601) Crossref, Web of Science, Google Scholar - 47.
Erlander S, Stewart NF . 1990 The gravity model in transportation analysis – theory and extensions. Zeist, The Netherlands: VSP Books. Google Scholar - 48.
Truscott J, Ferguson NM . 2012 Evaluating the adequacy of gravity models as a description of human mobility for epidemic modelling. PLoS Comput. Biol. 8, e1002699. (doi:10.1371/journal.pcbi.1002699) Crossref, PubMed, Web of Science, Google Scholar - 49.
Belton V, Stewart T . 2003 Multiple criteria decision analysis – an integrated approach. Dordrecht, The Netherlands: Kluwer Academic Publishers. Google Scholar - 50.
M’ikanatha N, Iskander J . 2015 Concepts and methods in infectious disease surveillance. Chichester, UK: Wiley Blackwell. Google Scholar - 51.
Barzilay E et al. 2013 Cholera surveillance during the Haiti epidemic – the first 2 years. New Engl. J. Med. 368, 599-609. (doi:10.1056/NEJMoa1204927) Crossref, PubMed, Web of Science, Google Scholar - 52.
Santa-Olalla P et al. 2013 Implementation of an alert and response system in Haiti during the early stage of the response to the cholera epidemic. Am. J. Trop. Med. Hyg. 89, 688-697. (doi:10.4269/ajtmh.13-0267) Crossref, PubMed, Web of Science, Google Scholar - 53.
Eppinger M et al. 2014 Genomic epidemiology of the Haitian cholera outbreak: a single introduction followed by rapid, extensive, and continued spread characterized the onset of the epidemic. mBio 5, e0172114. (doi:10.1128/mBio.01721-14) Crossref, Web of Science, Google Scholar - 54.
Frerichs R, Keim P, Barrais R, Piarroux R . 2012 Nepalese origin of cholera epidemic in Haiti. Clin. Microbiol. Infect. 18, E158-E163. (doi:10.1111/j.1469-0691.2012.03841.x) Crossref, PubMed, Web of Science, Google Scholar - 55.
Gaudart J, Rebaudet S, Barrais R, Boncy J, Faucher B, Piarroux M, Magloire R, Thimothe G, Piarroux R . 2013 Spatio-temporal dynamics of cholera during the first year of the epidemic in Haiti. PLoS Negl. Trop. Dis. 7, e2145. (doi:10.1371/journal.pntd.0002145) Crossref, PubMed, Web of Science, Google Scholar - 56.
Orata F, Keim P, Boucher Y . 2014 The 2010 cholera outbreak in Haiti: how science solved a controversy. PLoS Pathog. 10, e1003967. (doi:10.1371/journal.ppat.1003967) Crossref, PubMed, Web of Science, Google Scholar - 57.
Jutla A, Whitcombe E, Hasan N, Haley B, Akanda A, Huq A, Alam M, Sack RB, Colwell RR . 2013 Environmental factors influencing epidemic cholera. Am. J. Trop. Med. Hyg. 89, 597-607. (doi:10.4269/ajtmh.12-0721) Crossref, PubMed, Web of Science, Google Scholar - 58.
Righetto L, Casagrandi R, Bertuzzo E, Mari L, Gatto M, Rodriguez-Iturbe I, Rinaldo A . 2012 The role of aquatic reservoir fluctuations in long-term cholera patterns. Epidemics 4, 33-42. (doi:10.1016/j.epidem.2011.11.002) Crossref, PubMed, Web of Science, Google Scholar - 59.
Vesipa R, Ridolfi L . 2017 Impact of seasonal forcing on reactive ecological systems. J. Theor. Biol. 419, 23-35. (doi:10.1016/j.jtbi.2017.01.036) Crossref, PubMed, Web of Science, Google Scholar - 60.
Rinaldo A, Bertuzzo E, Blokesch M, Mari L, Gatto M . 2017 Modeling key drivers of cholera transmission dynamics provides new perspectives on parasitology. Trends Parasitol. 33, 587-599. (doi:10.1016/j.pt.2017.04.002) Crossref, PubMed, Web of Science, Google Scholar - 61.
Ciddio M, Mari L, Sokolow SH, De Leo G, Casagrandi R, Gatto M . 2017 The spatial spread of schistosomiasis: a multidimensional network model applied to Saint-Louis region, Senegal. Adv. Water Res. 108, 406-415. (doi:10.1016/j.advwatres.2016.10.012) Crossref, PubMed, Web of Science, Google Scholar - 62.
Mari L, Gatto M, Ciddio M, Dia ED, Sokolow SH, De Leo G, Casagrandi R . 2017 Big-data-driven modeling unveils country-wide drivers of endemic schistosomiasis. Sci. Rep. 7, 489. (doi:10.1038/s41598-017-00493-1) Crossref, PubMed, Web of Science, Google Scholar - 63.
Perez-Saez J et al. 2015 A theoretical analysis of the geography of schistosomiasis in Burkina Faso highlights the roles of human mobility and water resources development in disease transmission. PLoS Negl. Trop. Dis. 9, e0004127. (doi:10.1371/journal.pntd.0004127) Crossref, PubMed, Web of Science, Google Scholar - 64.
Xu B, Gong P, Seto E, Liang S, Yang C, Wen S, Qiu D, Gu X, Spear R . 2006 A spatial-temporal model for assessing the effects of intervillage connectivity in schistosomiasis transmission. Ann. Assoc. Am. Geogr. 96, 31-46. (doi:10.1111/j.1467-8306.2006.00497.x) Crossref, Google Scholar