Heterogeneity in schistosomiasis transmission dynamics

Highlights • Transmission dynamics of schistosomiasis presents multiple heterogeneity sources.• A comprehensive framework for heterogeneous disease transmission is proposed.• Heterogeneous multigroup communities can be more prone to parasite transmission.• Presence of multiple water sources can hinder parasite transmission.• Spatial and temporal heterogeneities can have nontrivial implications for endemicity.


Introduction
Transmission heterogeneity is a central issue in the study of infectious disease dynamics; in fact, it has been receiving considerable attention for more than thirty years, with research focusing on the role played by heterogeneity in both short-term epidemic dynamics and long-term transmission maintenance, as well as on the challenges and opportunities that heterogeneity poses to disease control (see e.g. Andreasen and Christiansen, 1989;Bolzoni et al., 2007;Diekmann et al., 1990;Dushoff and Levin, 1994;Heffernan et al., 2005;Hethcote and Van Ark, 1987;Lloyd-Smith et al., 2005;Nold, 1980;van den Driessche and Watmough, 2002;Woolhouse et al., 1997 ; see also VanderWaal and Ezenwa, 2016 for a recent review). A full understanding of the drivers and the consequences of heterogeneity still represents a major challenge for epi-demiology, especially for diseases characterized by indirect transmission, such as vector-borne, water-related, or macroparasitic infections ( Hollingsworth et al., 2015 ).
One such disease for which there is a unanimous scientific consensus about the relevance of heterogeneity on transmission is schistosomiasis, one of the most common parasitic diseases (second only to malaria) and the deadliest among the neglected tropical diseases ( Centers for Disease Control and Prevention, 2011 ). Schistosomiasis is caused by trematode parasites belonging to the genus Schistosoma , which need some species of freshwater snails as obligate intermediate hosts. It affects more than 250 million people in tropical and subtropical regions of the developing world, especially in sub-Saharan Africa, which is home to approximately 90% of worldwide cases ( World Health Organization, 2017 ). Human infection occurs through contact with water contaminated with freely swimming, short-lived schistosome larvae, known as cercariae, which are shed by infected snails. Specifically, cercariae can infect exposed humans through skin penetration. Within human hosts, they develop into sexually mature schistosomes. Depending on their species, the resting place of adult parasites is in the blood vessels around the bladder or the intestine of the infected human host. There, schistosomes mate and produce eggs that are excreted through urine or faeces. After reaching freshwater, eggs hatch into so-called miracidia, a second short-lived larval form of the parasite that can infect species-specific snail hosts and undergo asexual replication therein. Infected snails complete the parasite's life cycle by shedding cercariae into water (see e.g. Colley et al., 2014 ).
The transmission of schistosomiasis is thus controlled by contact with environmental freshwater infested with parasite larvae. Unsurprisingly, then, the disease is especially prevalent in communities that lack access to piped drinking water and adequate sanitation, and that have to resort to unsafe water sources for their primary needs ( Rollinson et al., 2013 ). Conversely, the availability of adequate water provisioning and sanitation infrastructures may offer protection against schistosomiasis ( Grimes et al., 2014 ). At a regional scale, different communities may be endowed with differential infection risk linked to their geographical and/or socioeconomic context. However, safe water supplies cannot completely avert human contact with environmental freshwater, nor can the presence of adequate sanitation guarantee per se its use ( Grimes et al., 2015 ). Agricultural, domestic, occupational and recreational tasks may foster contact with potentially infested water, and thus represent risk factors even in places where water provisioning and sanitation are adequate. Also, differences in behavior and lifestyle represent major potential determinants of transmission heterogeneity . As a result, the risk of schistosomiasis infection is usually not equally spread even at a local level, i.e. within a community, actually being higher for people whose routinary activities bring them in contact with water, such as fishermen, farmers and women. Overall, however, infection risk is highest among school-aged children, mostly because they are likely to spend a considerable amount of time in unsafe water; also, their immune system may not be fully developed yet ( World Health Organization, 2017 ). Several risk groups may thus be identified within a local community, possibly based on age, sex and/or occupation (e.g. Spear et al., 2004a ).
In communities with access to several water sources, even individuals within the same risk group may be endowed with differential infection risk because of their specific water-contact patterns. Different water bodies may represent habitats of different quality for the intermediate snail hosts of the schistosomes, and can thus be relatively heterogeneous in terms of the local snail population abundance they can host (e.g. Gurarie and King, 2005 ). As such, a comprehensive description of heterogeneity in schistosomiasis transmission also requires an account of individual preferences related to water contacts (e.g. Scott et al., 2003 ), as well as an ecological assessment of snail density (and possibly infection prevalence) at each of the available water sources ( Perez-Saez et al., 2016;Gurarie et al., 2017 ). Additionally, infection risk may also be heterogeneously distributed over time. Temporal variability of schistosomiasis transmission can arise from the seasonal fluctuations of climatological variables (e.g. Sturrock et al., 2001;Liang et al., 2007 ) -most notably rainfall/temperature in tropical/temperate areas ( McCreesh and Booth, 2013 ) -and environmental conditions (as in the case of fluctuating transmission risk during floods; e.g. Wu et al., 2008 ), both of which can remarkably influence snail population ecology (e.g. Chandiwana et al., 1987;Woolhouse et al., 1998;Sturrock et al., 2001 ). Water-contact patterns can also vary over time as a response to seasonal changes in human activities related, for instance, to agriculture and farming (e.g. Spear et al., 2004b;Liang et al., 2007 ). As a result of seasonal environmental forcing, transmission intensity and infection prevalence in both human and snail hosts are often observed to follow seasonal patterns (e.g. Chandiwana et al., 1987;Woolhouse et al., 1989;Sturrock et al., 2001 ).
Heterogeneity has been accounted for in mathematical models for schistosomiasis transmission since long. A seminal contribution in this respect was given by Barbour (1978) , who elaborated on Macdonald 's (1965) basic model, and extended it to account for individual variations in water contact patterns and source heterogeneity. In particular, Barbour derived an expression for the threshold parameter controlling disease transmission and concluded that heterogeneities in water contacts can favor endemicity. Later works focused on the role played by multiple water sources and individual heterogeneity ( Woolhouse et al., 1998, site-specific environmental features ( Liang et al.,20 02,20 05 ), and the interplay between local-scale heterogeneities and spatial coupling mechanisms ( Gurarie and King, 2005 ) in schistosomiasis transmission dynamics, with a special emphasis for the integration of field data and the analysis of strategies for disease control. Recent modeling studies have also focused on spatial heterogeneity arising from the pathogens' dispersal mechanisms (such as human mobility for adult schistosomes and hydrological transport for their larval stages), as well as on the implications for parasite endemism (e.g. Gurarie and Seto, 2009;Perez-Saez et al., 2015;Ciddio et al., 2017 ). Although less frequently than transmission/spatial heterogeneity, temporal heterogeneity in schistosomiasis transmission has sometimes been also accounted for in both applied (e.g. Liang et al., 20 02, 20 05;Remais, 2010;Gurarie et al., 2017 ) and theoretical ( Zhang et al., 2014;Ciddio et al., 2015 ) studies.
In this work we aim to provide a comprehensive -yet assimple-as-possible -account of the major sources of heterogeneity that can be relevant to schistosomiasis transmission. To that end, we make use of suitable extensions of the seminal Macdonald 's (1965) model, whose main characteristics are briefly reviewed in the next section. Then, we analyze a general transmission model (after Barbour, 1978 ) accounting for heterogeneities in both the human host population (expressed as differential exposure/contamination risk associated with different sub-groups within a community) and in the available water sources (as resulting from specific water contact patterns and the distribution of the snail host population). A threshold condition for the multigroup/multi-source model is worked out based on the dominant eigenvalue of a generalized reproduction matrix accounting for multiple sources of heterogeneity. The cases of heterogeneity in the human host population and in the available water sources are then analyzed separately through simple numerical examples, so as to single out the specific effects of the different sources of heterogeneity on schistosomiasis transmission dynamics. The case of spatial heterogeneity (which can be thought of as a special combination of mixed host/source heterogeneity) is also analyzed in some detail. Finally, a seasonally-forced version of the basic Macdonald 's (1965) model is studied, and conditions for long-term parasite establishment in periodic environments are obtained. A discussion about the most important epidemiological consequences of heterogeneity in schistosomiasis dynamics, with special focus on application to field data and disease control, closes the paper.

The basic homogeneous model
The simplest model for schistosomiasis was proposed by Macdonald in 1965 . It describes transmission dynamics through two state variables, namely the average parasite burden in the human population ( P ) and the prevalence of infection in snails ( Y ). Population dynamics are neglected in human and snail hosts, assuming instead demographic equilibrium for both. Also, the model does not include the dynamics of cercariae and miracidia, whose abundances are considered to be proportional to Y and P , respectively. The model thus reads as follows: where β is the snail-to-human transmission rate, N is the abundance of the snail population, γ is the mortality rate of adult parasites in human hosts, χ is the human-to-snail transmission rate, H is the abundance of the human population and μ is rate at which infected snails die and are replaced by uninfected ones. The parameters γ and μ can be evaluated as the inverse of the average lifespans of adult worms and infected snails (around five years and two months, respectively; see e.g. Feng et al., 2004 ), whereas β and χ represent aggregated parameters accounting for several epidemiological and socioeconomical processes.
Model (1) has two steady-state solutions. The first one is the so-called disease-free equilibrium (DFE), i.e. a state of the system in which the parasite is not present (thus P = 0 and Ȳ = 0 ). The second one is the endemic equilibrium (EE), i.e a state of the system in which parasite transmission is permanent, namely By standard linear stability analysis arguments (Appendix A online), it can be shown that the DFE is stable if γ μ − βNχH > 0 , also corresponding to the parameter region in which the EE is not feasible (i.e. characterized by negative components). This stability condition can be equivalently stated in terms of the so-called basic reproduction number i.e. the average number of established and reproductively mature offspring produced by a mature parasite during its lifetime in a population of uninfected hosts ( Anderson and May, 1992 ). Specifically, if r 0 < 1 the DFE is asymptotically stable, while the EE is unfeasible and unstable. Conversely, if r 0 > 1 the DFE is unstable, while the EE is feasible and stable (see again Appendix A). Fig. 1 illustrates the basic properties of the equilibria of model (1) . In particular, the shape of the contour lines for the average parasite burden in human hosts explains why preventing water contamination (i.e. decreasing the overall human-to-snail transmission rate, χH ) might be less an effective measure for disease control (namely for decreasing P ) than preventing human water contact (i.e. decreasing the overall snail-to-human transmission rate, βN ), as already suggested by Macdonald (1965) in his seminal work (see also Grimes et al., 2015 , for discussion).

A general model for heterogeneous schistosomiasis transmission
The basic Macdonald 's (1965) transmission model (1) does not account for any of the different sources of heterogeneity that may be relevant to schistosomiasis transmission. Specifically, it does not include heterogeneities arising, for instance, from the presence of sub-groups in the human population, each of which may be endowed with different exposure/contamination rates, and/or from the availability of several water sources, whose characteristics (inclusive of local snail abundances) may possibly yield differential transmission risk. A general, multi-group model for schistosomiasis transmission was proposed by Barbour in 1978 . Following his approach, we consider a community subdivided in G groups (e.g. according to age, sex or occupation, or to any other trait that may be relevant to schistosomiasis transmission) with access to S different water sources, all of which can possibly host snail populations. The model can be stated as Macdonald 's (1965) homogeneous model (1) . The DFE is stable if r 0 < 1 (gray-shaded parameter combinations), while the EE is feasible and stable if r 0 > 1. The black solid line indicates the transcritical bifurcation through which the two equilibria collide and exchange stability (endemicity boundary, r 0 = 1 ). The gray dashed curves are contour lines for the prevalence of infected snails, while the gray dash-dotted curves refer to the average worm burden in human hosts. Parameter values: γ = 5 . 5 · 10 −4 , μ = 1 . 7 · 10 −2 . All rates are expressed in [day −1 ].

dY s dt
where the subscript g [ s ] indicates variables or parameters pertaining to different human groups [water sources], while E gs [ C gs ] represents one entry of the human exposure [contamination] matrix, which describes the specific snail-to-human [human-to-snail] transmission risk associated with people from group g being in contact with water at source s . Clearly, a fully mechanistic description of the interwined processes that ultimately determine actual exposure/contamination risk for humans may be -at best -impractical. For the sake of simplicity, we thus describe human exposure/contamination as the product of three terms each, namely ] is a group-specific relative exposure [contamination] risk, ω gs represents the fraction of water contacts occurring at source s for people belonging to group g (so that 0 ≤ ω gs ≤ 1 and s ω gs = 1 ) and φ E s [ φ C s ] is a source-specific relative factor contributing to exposure [contamination]. Group-specific scores reflect differences in transmission risk related to demographic and socioeconomic traits. As an example, groups consisting prevalently of fishermen vs. blue-or white-collar workers may be endowed with different values of E g and C g . Source-specific scores, instead, quantify differential transmission risk associated with the environmental characteristic of the water bodies. For instance, a shallow ephemeral pond may have intrinsically different values of φ E s and φ C s than a permanent river site. From a technical perspective, we assume that the union of the graphs associated with the exposure and contamination matrices ( E = [ E gs ] and C = [ C gs ] , respectively) is strongly connected, so as to prevent human groups and/or water sources from possibly being isolated as far as water contact is concerned. Equilibrium solutions of model (2) include a DFE, in which P g = 0 for all groups and Ȳ s = 0 for all sources, and an EE, whose components can be determined analytically only in some particular cases ( Barbour, 1978 ). It can be shown (Appendix B) that the stability of the DFE of model (2) depends on the generalized reproduction matrix in which h and n are diagonal matrices whose non-zero elements correspond, respectively, to the fraction h g = H g /H of the total human host population ( H ) belonging to each sub-group ( G g=1 h g = 1 ) and to the fraction n s = N s /N of the overall snail population ( N ) that can be found at each water site ( S s =1 n s = 1 ), and the superscript T indicates matrix transposition. Specifically, the DFE of model (2) is asymptotically unstable if the dominant eigenvalue of matrix R ( R GS cates that the threshold parameter refers to a system with G human groups and S water sources) is larger than one, stable otherwise. This condition thus generalizes the analogous threshold based on the basic reproduction number r 0 that is found for the basic Macdonald (1965) model (see Barbour, 1978 ;Gurarie and King, 2005) . If R GS 0 > 1 , numerical simulations of model (2) show that the EE is stable and characterized by strictly positive components.
In general, establishing parasite invasion/establishment conditions in the presence of transmission heterogeneity requires the numerical evaluation of R GS 0 . In the next sections we analyze two important special cases of model (2) for which the threshold parameter can be derived analytically.

Group heterogeneity: G human groups, one water source
We first analyze the case in which G ( > 1) sub-groups can be identified in a human population with access to one single water source ( S = 1 ). In this case, ω = [ ω gs ] = i G , where i G is a column vector of length G with all elements equal to one. With no loss of generality, we can assume φ E s = φ C s = 1 , as changes in the source-specific exposure/contamination risk factors can be absorbed -with one water source -into the transmission rates β and χ. To allow a fair comparison with the results of the homogeneous model (1) , we assume that the total human population abundance evaluated over all groups and the average exposure/contamination risk of the heterogeneous multi-group community are the same of an equivalent homogeneous community (i.e. g h g = 1 and g h g E g = g h g C g = 1 ). In the presence of group heterogeneity, the DFE is unstable, thus allowing for parasite invasion and long-term establishment in the community, if As already noted by Woolhouse et al. (1998) , under the additional hypothesis that for each group the relative exposure risk is equal to the relative contamination risk ( E g = C g = g ), it is possible to show that heterogeneous multi-group communities sharing one water source cannot be less prone to long-term parasite establishment than homogeneous ones, i.e. that R G 1 0 ≥ r 0 , all else being equal (see again Appendix C). Indeed, R G 1 0 = r 0 only in the case of a homogeneous multi-group community ( g = 1 for all groups). Note that the simplification E g = C g may indeed be reasonable for a disease, like schistosomiasis, in which exposure and contamination are both tightly related to water contact Woolhouse et al., 1998;Gurarie and Seto, 2009 ). People spending more time at water sources where they can be exposed to the pathogen are also more likely to contribute to water contamination with their excreta. Therefore, the same behavior that increases exposure risk is also likely to increase contamination risk.
If R G 1 0 > 1 , the DFE is unstable and endemic schistosomiasis transmission is possible (stable EE). Three noteworthy results emerge for the EE of a heterogeneous multi-group community with one water source (details in Appendix C): a) the prevalence of infected snails cannot be smaller than in an equivalent homogeneous community ( Ȳ G 1 ≥Ȳ ); b) even some groups whose exposure/contamination risk is not larger than that of an equivalent homogeneous community ( g ≤ 1) may have a parasite load not smaller than that of an equivalent homogeneous community ( P G 1 g ≥P for g ≥ , with ≤ 1); and c) the average parasite burden is not smaller than in an equivalent homogeneous community ( g h g P G 1 g ≥P ). These three results depend on the fact that g h g 2 g ≥ 1 in a heterogeneous multi-group community with one water source.

The case of two groups
To elucidate the implications of heterogeneity associated with differential transmission risk, we analyze the simplest case of between-group heterogeneity, namely a community in which human hosts can be divided in two sub-populations (i.e. G = 2 ), and where for each group the relative exposure risk is equal to the relative contamination risk ( E g = C g = g ). Let f be the fraction of human hosts belonging to the higher-risk group (say group 1, h 1 = f and h 2 = 1 − f, 0 < f < 1) and let k be the relative exposure/contamination risk of group 1 compared to group 2 ( 1 / 2 = k, k ≥ 1). To bettere contrast the effects of heterogeneous vs. homogeneous transmission, it is convenient to impose that the average exposure/contamination risk in the heterogeneous two-group community is the same as in the homogeneous one, that is h 1 1 + h 2 2 = 1 , thus obtaining A few quantitative examples can help elucidate to what extent a two-group heterogeneous community can be more prone to the establishment of endemic schistosomiasis transmission than an equivalent homogeneous community. For instance, if a small fraction of the population (say 10%) is characterized by relatively large exposure/contamination risk (say 10 times higher) compared to the rest of the population, then R 21 0 ≈ 3 r 0 . If heterogeneity is even more pronounced (say 0.5% of the population with a 500fold higher infection risk), then the parasite can establish in the whole community (i.e. also in the lower-risk group, albeit with a low average burden) even if r 0 is as low as ≈ 0.01. Fig. 2 A generalizes these results and shows that a common feature of heterogeneous multi-group communities is sub-threshold endemic transmission ( van den Driessche and Watmough, 2002 ), i.e. the longterm establishment of the parasite ( R 21 0 > 1 ) under conditions that would prevent it in an equivalent homogeneous community ( r 0 < 1, for which P = 0 and Ȳ = 0 ). Panel B of Fig. 2 shows the steadystate components of the EE in a heterogeneous two-group community in the case of sub-threshold transmission ( r 0 = 0 . 9 ). The three theoretical results described above (namely, Ȳ 21 ≥Ȳ , P 21 2 ≥P for 2 ≤ 1, and g h g P 21 g ≥P ) are immediately recovered. Furthermore, we note that both the average parasite burden in the highrisk group and the prevalence of infected snails can attain very large values in case of pronounced heterogeneity. Not surprisingly, the average parasite burden in the low-risk group can be much (indeed, k times) smaller than in the high-risk one; however, because of the different group sizes, most parasites may actually be hosted within the members of the low-risk group.
More complex scenarios, i.e. describing communities with several (possibly many, G 2) risk groups relative to exposure and/or  (2) with G = 2 and S = 1 ). The two groups differ for their relative abundance ( h 1 and h 2 , with h 1 + h 2 = 1 ) and intrinsic transmission risk ( 1 and 2 ). A) Endemicity boundaries: parasite establishment is possible ( R 21 0 > 1 ) for parameter combinations lying above the bifurcation curves, which correspond to R 21 0 = 1 and are obtained for different values of the basic reproduction number r 0 of an equivalent homogeneous community (labels). B) State variables at the stable equilibrium: the dashed curves are contour lines for the prevalence of infected snails, while the dash-dotted curves refer to the average worm burden in human hosts (dark gray: high-risk group; light gray: low-risk group). The DFE is stable if R 21 0 < 1 (gray-shaded parameter combinations), while the EE is feasible and stable if R 21 0 > 1 . The black solid line indicates the endemicity boundary R 21 0 = 1 . In panel B, the overall exposure ( βN ) and contamination ( χ H ) rates are assumed to be equal, and their value is set so as to match that of a homogeneous community endowed with a basic reproduction number r 0 = 0 . 9 (thus, βN = χ H = √ r 0 γ μ). Other parameters as in Fig. 1 .
contamination risk and a single water source are discussed in Appendix C.

Source heterogeneity: homogeneous human population, S water sources
We now turn to the study of a homogeneous human community ( G = 1 ) with access to a set of S ( > 1) different water sources.
In this case, matrix ω reduces to a row vector u = [ ω 1 · · · ω S ] , whose elements ω s represent the fraction of water contacts that take place at each source s . With no loss of generality, we can set E g = C g = 1 , as changes in group-specific exposure/contamination risk factors could be incorporated -with just one group -into the transmission rates β and χ. For the sake of parameter parsimony, and to avoid possible confounding effects, we also assume φ E s = φ C s = 1 , i.e. that water sources share the same characteristics leading to exposure/contamination, including their accessible water volume. On the other hand, heterogeneity is granted by the distribution of the overall snail population among different water sources, as well as by the distribution of water contacts. To allow a fair comparison with the homogeneous case of model (1) , the total snail abundance and the overall frequency of human-water contacts evaluated over all water points are assumed to be the same as in an equivalent homogeneous system (i.e. s n s = s ω s = 1 ).
In the presence of source heterogeneity, the DFE is unstable and endemic transmission is possible if Details are reported in Appendix D, where it is also shown that R 1 S 0 ≤ r 0 for all water contact and snail distribution patterns. Therefore, as already noted by Woolhouse et al. (1991) , a dilution effect induced by the parcelling of the overall snail population and the human water contacts among the available sources exists. This dilution effect can decrease the risk of schistosomiasis endemicity, all else being equal. If R 1 S 0 > 1 , the DFE is unstable. In this case, although a compact analytical expression for its coordinates is not readily available, the stability of the EE can easily be tested via numerical simulation.

The case of two water sources
To better analyze parasite establishment in a community with access to S water sources we refer to the simplest setting, namely the case in which S = 2 . Fig. 3 illustrates the dilution effect induced by the presence of multiple water sources, in terms of both infection intensity and endemicity risk. Specifically, panels A-C shows a comparison between the equilibrium values of the state variables describing the heterogeneous two-source community at equilibrium (obtained via numerical simulation) and the equivalent homogeneous case of model (1) . Both the prevalence of infection in snails and the average parasite burden in human hosts are found to be smaller in the two-source community. In particular, superthreshold parasite extinction (i.e. occurring for r 0 > 1) is possible if a large share of water contacts occur at a water source hosting a small fraction of the overall snail population (or equivalently, if a small share of water contacts occur at a water source hosting a large fraction of the overall snail population; gray-shaded parameter combinations). Panel D of Fig. 3 shows instead that R 1 S 0 approaches r 0 only if n 1 and ω 1 are either both large or both small, i.e. in cases in which heterogeneity is practically negligible. Conversely, R 1 S 0 progressively decreases compared to r 0 as more water contacts are made at relatively snail-free sources (i.e. if n 1 is large and ω 1 is small or, viceversa, if n 1 is small and ω 1 is large). In case of uniform snail distribution ( n 1 = n 2 = 1 / 2 ), endemicity risk decreases (higher values of r 0 are needed for the DFE to become unstable) as the heterogeneity of human-water contact also decreases ( | ω 1 − ω 2 | → 0 ). Similar analyses can be performed with reference to more complex case studies, i.e. communities with access to several (possibly many, S 1) water sources. Some examples are discussed in Appendix D.

Spatial heterogeneity: multiple goups and water sources
Model (2) can be used to study the effects of heterogeneity also in a spatially-explicit setting. The G human groups could in fact belong to different villages; also, they could have access to a set of S water sources arranged in a network endowed with a well-defined spatial structure (see e.g. Gurarie   Analysis of a single-group community with access to two water sources (model (2) with G = 1 and S = 2 ). The two sources differ for the relative abundance of snails they host ( n 1 , bottom axes, and n 2 , top axes, with n 1 + n 2 = 1 ) and the frequency of human-water contacts ( ω 1 , left axes, and ω 2 , right axes, with ω 1 + ω 2 = 1 ). A) Steady-state prevalence of infection in snails at water source 1 ( Ȳ 12 1 ) divided by the steady-state prevalence of infected snails in an equivalent homogeneous community ( Ȳ ) with r 0 = 3 ( βN = χ H = √ r 0 γ μ, other parameters as in Fig. 1 ). Black curves are contour lines of Ȳ 12 1 / Ȳ . B) As in panel A, for the steady-state prevalence of infection in snails at water source 2 ( Ȳ 12 2 ). Black curves are contour lines of Ȳ 12 2 / Ȳ . C) As in panel A, for the steady-state average parasite burden in human hosts ( P 12 ). Black curves are contour lines of P 12 / P (where P is the steady-state average parasite burden in the equivalent homogeneous community). The 0-level lines in panels A-C (thick lines) correspond to the endemicity boundary R 12 0 = 1 , while the gray-shaded regions correspond to parameter combinations for which the DFE of the heterogeneous model is stable ( R 12 0 < 1 ). D) Endemicity boundaries in the heterogeneous community ( R 12 0 = 1 ) for different values of the basic reproduction number of an equivalent homogeneous community (labels). Parasite establishment is possible ( R 12 0 > 1 ) for parameter combinations lying within the region(s) including the homogeneous cases n 1 = ω 1 = 0 and/or n 1 = ω 1 = 1 .
water contact matrix ω can be thought of as a human mobility matrix, at least as far as water use is concerned. Here we limit our attention to the simplest case of spatial heterogeneity, namely a metacommunity with two sub-groups ( G = 2 ) living in two geographically distinct locations, each of which has one preferential water source ( S = 2 ). If we assume that water contacts occur not only at the source that is closest to the home site, but also when people travel between the two sites, the water contact matrix can be defined as where the anti-diagonal elements represent the fraction of water contacts occurring away from the home site. The only spatial coupling process considered in the model is human mobility, which seems to be appropriate if water sources are not hydrologically connected, or if they are quite distant from each other.
To allow a fair comparison with the (spatially) homogeneous case, we use the same human and snail total population abun-dances as in the homogeneous case. Furthermore, we assume that E i = C i = i and φ E i = φ C i = 1 for both i = 1 , 2 , and define 1 and 2 so that the average exposure/contamination risk in the (spatially) heterogeneous metacommunity is the same as in the homogeneous community. Finally, we set m 1 = m 2 = m for the sake of simplicity. The stability of the DFE ([0 0 0 0] T ) can be assessed through the evaluation of the dominant eigenvalue of matrix R , which can be easily performed numerically. However, in spite of its simplicity, the spatially explicit version of the transmission model still holds many degrees of freedom, namely: the values of the average exposure and contamination rates (whose product concurs to the definition of the basic reproduction number in the equivalent homogeneous community); the level of human mobility, as measured by the fraction of water contacts occurring away from the home site; the distributions of the human/snail populations in the two sites; and the spatial heterogeneity of exposure/contamination risk. This makes a complete study of the endemicity threshold quite cumbersome. Fig. 4 thus reports the analysis of parasite establishment conditions for four selected settings. Fig. 4. Endemicity boundaries in a spatially heterogeneous metacommunity with two human groups living in separate locations, each with one preferential water source (model (2) with G = 2 and S = 2 ). The groups may differ for their relative abundance ( h 1 and h 2 , with h 1 + h 2 = 1 ) and intrinsic transmission risk ( 1 and 2 ), while the two sources may differ for the relative abundance of snails they host ( n 1 and n 2 , with n 1 + n 2 = 1 ) and the frequency of human-water contacts (for each community, the fraction of contacts at the farthest water point is m , while 1 − m is the fraction of contacts at the home site). Parasite establishment is possible ( R 22 0 > 1 ) on the right of the bifurcation curves, which correspond to R 22 0 = 1 and are obtained for different spatial distributions of the human host population (legend). A) Spatially homogeneous transmission risk ( 1 / 2 = 1 ) and snail population distribution ( n 1 = n 2 = 1 / 2 ). B) As in panel A, with spatially heterogeneous transmission risk ( 1 / 2 = 10 ). C) As in panel A, with spatially heterogeneous snail abundance ( n 1 = h 1 ). D) Spatially heterogeneous transmission risk ( 1 / 2 = 10 ) and snail population distribution ( n 1 = h 1 ).
In the first example ( Fig. 4 A), exposure/contamination risk is the same in the two sites and the snail population is uniformly distributed between the two water sources. The endemicity boundary is symmetric with respect to h 1 = h 2 = 1 / 2 and m = 1 / 2 , which means that the values of R 22 0 evaluated for h i = h or h i = 1 −h , and for m = m or m = 1 −m coincide. Therefore, only results obtained for 0 ≤ h 1 ≤ 1/2 and 0 ≤ m ≤ 1/2 are shown. In this case, mobility shows a well-defined dilution effect, with increasing m values playing against long-term pathogen establishment. The endemicity boundary generally (i.e. except for m = 1 / 2 , corresponding to complete spatial mixing) moves rightward as the fraction of residents of site 1 increases from 0 to 1/2, which indicates that more homogeneous population distributions are less conducive to the establishment of endemic parasite transmission. For instance, if the total human population were uniformly distributed between the two sites, the parameters combinations yielding R 22 0 > 1 would give r 0 > 4 (note that, in this case, each village/water source would host exactly half of the human/snail population of an equivalent homogeneous community).
The second example ( Fig. 4 B) accounts for spatial heterogeneity in terms of exposure/contamination risk. This might represent the (quite common) case of a village composed by different neighborhoods that share the same water contact points but are composed of a diversified human population, with people characterized by distinct socioeconomic traits. Such differences in demographic composition could indeed yield a correspondence between spatial and transmission heterogeneities, as mediated by different water-related behavior (which would most likely be preserved during movement, especially over relatively short spatial scales). For instance, Pinot de Moira et al. (2007) found an approximately 10fold difference in the overall duration of water contacts (mostly related to fishing activities) between groups of adult males with different tribal backgrounds living in separated neighborhoods of a community in rural Uganda. Because of the spatial heterogeneity in exposure/contamination, the endemicity boundary is no longer symmetric about h 1 = h 2 = 1 / 2 , while it still is about m = 1 / 2 . Therefore, only results obtained for 0 ≤ m ≤ 1/2 are shown. Interestingly, it turns out that R 22 0 is maximum compared to r 0 for intermediate values of h 1 (actually, for 0 < h 1 < 0.2). In other words, the likelihood of endemicity is maximum if a relatively small fraction ( ≈ 10% in this example) of the human population lives in the village where the most-at-risk water source is located. Mobility is again associated with a dilution effect. For some population distributions and mobility levels (i.e. for h 1 = 0 . 1 and m = 0 . 05 ), R 22 0 can indeed be larger than one even if r 0 is not (sub-threshold parasite establishment), which suggests that the spatial heterogeneity of the transmission parameters can strongly favor endemic transmission also in the presence of human mobility.
The third example ( Fig. 4 C) differs from the first one in that snail population abundance is spatially heterogeneous. Specifically, the spatial distribution of snail hosts is assumed to be linked to that of the human population ( n i = h i ). This might be reasonable, for instance, should a positive feedback between larger human groups and larger snail populations exist, as it may be expected in rural contexts where man-made alterations of the environment (such as the construction of reservoirs and irrigation schemes for agricultural development) can contribute to creating new habitat for the snails ( Steinmann et al., 2006 ). Because of the spatial heterogeneity of snail distribution, the endemicity boundary is no longer symmetric with respect to m = 1 / 2 , while it still is about h 1 = h 2 = 1 / 2 . Therefore, only results obtained for 0 ≤ h 1 ≤ 1/2 are shown. High levels of human mobility are typically associated with parasite extinction (possibly except for h 1 = h 2 = 1 / 2 ). It is also interesting to note that the endemicity boundary moves rightward for increasing values of h 1 = n 1 only if m < 1/2, leftward otherwise ( m > 1/2).
In the fourth example ( Fig. 4 D) both the exposure/contamination risk and the snail distribution are assumed to be spatially heterogeneous. As a result, the endemicity boundary is symmetric about neither h 1 = h 2 = 1 / 2 nor m = 1 / 2 . The most interesting result pertaining to this case is the appearance of backward-bending endemicity boundaries. Of particular importance is the region identified by m < 1/2, where mobility is found to play a significant role. Specifically, for 0 < h 1 < 1/2, R 22 0 is minimum for intermediate values of m . Therefore, depending on h i and r 0 , increasing human mobility can either hinder or promote endemic parasite transmission. With h 1 = 0 . 2 and r 0 = 2 , for example, parasite establishment is possible for m = 0 or m = 0 . 4 , but not for m = 0 . 2 . Interestingly, even relatively low values of the 1 / 2 ratio are sufficient to recover qualitatively similar results (e.g. 1 / 2 = 2 ; see Appendix E for some examples). Finally, we note that cases in which more abundant snail populations are associated with less abundant human communities (e.g. n i = 1 − h i ) might instead be found in different socioeconomic contexts and/or spatial scales, such as in the case of a spatially explicit setting encompassing both urban and rural environments. Note that results for this latter scenario can be recovered from the former's, as the values of R 22 0 evaluated for m = m under the assumption n i = 1 − h i would coincide with those evaluate for m = 1 −m with n i = h i .

A special case of heterogeneity: time-varying transmission
So far we have been dealing with heterogeneity arising from two basic factors, namely differential infection risk among different sub-groups within a human community, and/or differential transmission related to the distribution of water contacts and snail hosts among the available water sources -with spatial heterogeneity representing a particular (yet noteworthy) combination of these two factors. In all cases, we have (mostly) focused on the relationship between transmission heterogeneity and parasite invasion/establishment. One possible source of heterogeneity that is not accounted for in models (1) and (2) (hence also in the special cases analyzed above) is the temporal variability of the parameters that are relevant to schistosomiasis dynamics. Some seasonallyforced models for schistosomiasis transmission have recently been proposed. For instance, Zhang et al. (2014) analyzed a model in which the human population is split into susceptible/infected hosts, without accounting for the average parasite burden in the community, as it is instead customary in models for macroparasite transmission. For this reason, their approach lies somewhat outside Macdonald 's (1965) modeling framework. Ciddio et al. (2015) , instead, proposed a Macdonald-like model accounting for the demographic dynamics of the snail population. They showed that the interplay between density dependence and temporal forcing may lead to complex (and realistic) transmission dynamics, yet they did not specifically focus on the conditions leading to endemic schistosomiasis transmission. Finally, other temporally-forced models (e.g. Liang et al., 20 02, 20 05;Remais, 2010;Gurarie et al., 2017 ) were mainly tailored to analyze real case studies, hence they are not fully amenable to theoretical investigation.
The effects of time-varying transmission on parasite invasion/establishment can be analyzed with the following timevarying version of model (1) , i.e.
To simplify the analysis of model (3) , we assume that the fluctuations of the transmission rates are mostly associated with seasonal variations of the human-water contact rate, so that ˜ β(t) and ˜ χ (t) are characterized by the same amplitude ( α β = α χ = α ω ) and phase ( ψ β = ψ χ = ψ ω ). Also, because of the long lifespan of schistosomes within human hosts, it is reasonable to neglect fluctuations of the mortality rate of adult parasites ( α γ = 0 ). Seasonal variations of the mortality rate of infected snails are not considered either ( α μ = 0 ), as they are better addressed in schistosomiasis transmission models that include a detailed description of the snail host ecology (see e.g. Ciddio et al., 2015;Sokolow et al., 2015;Gurarie et al., 2017 ). The period of oscillations is assumed to be one year for all the seasonal environmental signals ( τ x = τ = 365 [d], x ∈ { β, χ, N, H }). Without loss of generality, we can also set ψ ω = 0 and study the effects of possible lags ψ x between the phase of the transmission rates and the other parameters' (0 ≤ ψ N , ψ H ≤ τ ). A simplified version of model (3) incorporating all these assumptions is given in Appendix F.  (1) , the DFE of model (3) is a state in which P = 0 and Ȳ = 0 . However, determining the stability properties of the DFE (or, equivalently, the conditions under which the parasite can establish in the community) is relatively more complex in this case, because Floquet theory (instead of basic linear stability analysis) has to be applied (see e.g. Bacaër and Guernaoui, 2006 ;Klausmeier, 2008 ;Bittanti and Colaneri, 2009;Mari et al., 2014 ). Specifically, the DFE is unstable (thus allowing parasite invasion and the establishment of periodic tramsmission cycles) if and only if R F 0 , the maximum Floquet multiplier of system (3) linearized in a neighborhood of the DFE, is larger than one. Floquet multipliers are defined as the eigenvalues of the monodromy matrix M ( τ ), which can be obtained by solving the matrix differential equation

Like in model
over one period (i.e. over 0 ≤ t ≤ τ ), with the identity matrix as initial condition (see again Klausmeier, 2008 ). Note that R F 0 is equal to r 0 in the time-constant case ( α x = 0 for all x 's).
The results concerning parasite invasion in seasonal environments are reported in Fig. 5 . Seasonal fluctuations can make the instability of the DFE either less or more likely, depending on the parameter (panel A) or combination of parameters (panels B and C) that is assumed to vary over time. Specifically, if the abundance of one of the two host populations fluctuates periodically (0 < α N ≤ 1 or 0 < α H ≤ 1, e.g. because of seasonal variations of snail demography, or human migrations linked to seasonal activities), then R F 0 < r 0 . In these two cases, pathogen persistence requires larger average values of the time-varying parameters than in the time-constant case. On the contrary, R F 0 > r 0 if the human exposure/contamination varies seasonally (0 < α ω ≤ 1). In this latter case, endemic transmission is favored by seasonal forcing ( Fig. 5 A).
Coupled fluctuations of the human-water contact rate and the population abundance of either snails or humans are associated with the largest deviations from the endemicity threshold r 0 = 1 of the time-homogeneous case. Phase shifts between the two periodic signals are of paramount importance to qualify such deviations: on the one hand, synchronous fluctuations ( ψ N = 0 or ψ H = 0 , i.e. more frequent human-water contact when snails or human hosts are more abundant) promote sub-threshold endemicity, i.e. para-site invasion and establishment being possible for r 0 < 1; on the other hand, periodic signals that are in antiphase with each other ( ψ N = τ / 2 or ψ H = τ / 2 , e.g. more frequent human-water contact when snails or human hosts are less abundant) hinder endemicity, with parasite extinction being possible for r 0 > 1. Similar results are obtained with coupled fluctuations of the two hosts' population abundances, although with smaller deviations from the timehomogeneous case ( Fig. 5 B). The relative amplitude of fluctuations also matters in case of coupled seasonal forcing ( Fig. 5 C).

Discussion and conclusions
In this work we have analyzed a flexible modeling framework to describe schistosomiasis transmission dynamics in the presence of various sources of heterogeneity, including differential transmission risk linked to the demographic or socioeconomic traits of the human host population, availability of multiple water sources characterized by non-homogeneous water contacts and snail host distribution, spatially explicit water contact patterns, and temporal fluctuations of the eco-epidemiological parameters relevant to disease transmission. Overall, a clear picture emerges where heterogeneity plays a highly non-trivial role in the definition of epidemiological patterns. Specifically, between-group heterogeneity in transmission risk is typically associated with larger infection prevalence in snails and average parasite burden in humans, as well as with a higher likelihood of long-term parasite invasion. On the other hand, heterogeneity induced by the presence of different water sources may produce opposite effects, namely lower infection intensity in humans and snails, and lower chances of endemic parasite transmission compared to a homogeneous situation. All these results conform with and extend previous findings on heterogeneous schistosomiasis dynamics ( Barbour, 1978;Woolhouse et al., 1991Woolhouse et al., , 1998. Results concerning the spatially heterogeneous case reveal that the epidemiological implications of human mobility may be different in different contexts. Interestingly, there are cases in which the likelihood of endemic transmission is maximum for intermediate levels of human mobility. A similar scenario has been discussed in a recent modeling study of the geography of schistosomiasis transmission in Burkina Faso ( Perez-Saez et al., 2015 ). Therefore, despite the highly simplified setting used in this work, our analysis suggests that the interplay between transmission heterogeneity and spatial heterogeneity can provide a theoretical background to results found in real or realistic case studies. Fig. 6. Analysis of exemplificative strategies for schistosomiasis control in the presence of heterogeneity. Top panels refer to the case of a two-group community with access to a single water source (as in Fig. 2 with f = 0 . 1 and k = 10 ). A) Reduction of the reproduction number R 21 0 for increasing shares of the human population being involved in the intervention. B) Reduction of the average parasite burden in the two sub-groups (left axis: high-risk group; right axis: low-risk group). C) Reduction of infected snail prevalence. Bottom panels refer to the case of a single community with access to two water sources (as in Fig. 3 , with ω 1 = 3 / 4 , ω 2 = 1 / 4 and n 1 = n 2 = 1 / 2 ). D) Reduction of the reproduction number R 12 0 for increasing shares of the snail population being removed. E) Reduction of the average parasite burden. F) Reduction of infected snail prevalence at the two water sources. In all cases, βN = χ H are set so that the reproduction number in the heterogeneous community is equal to 2.5 prior to the intervention. Other parameters as in Fig. 1 . See text for details on the interventions.
Finally, the analysis of temporal heterogeneity in schistosomiasis transmission has revealed that the effects of seasonal forcing are context-dependent too. In fact, temporal heterogeneity can either promote or play against long-term parasite persistence in the community, depending on the eco-epidemiological parameters that are affected by environmental fluctuations. Although these results are noteworthy per se , the greater goals of mathematical epidemiology lie in the applicability of theoretical tools to real case studies and in the pursuit of opportunities for disease control. As far as parasite control is concerned, heterogeneity associated with differential transmission risk (group heterogeneity) seems to be the most critical and targetable factor that unambiguously favors schistosomiasis transmission. In this case, in fact, the presence of sub-threshold endemic dynamics suggests that it may be very difficult to break transmission in highly heterogeneous human communities, unless suitable interventions aimed at the most-at-risk share of the population (e.g. children, fishermen, or women, in the case of schistosomiasis) are implemented. The prototypical case of a two-group community with access to one single water source can be used again to get a handle on the importance of tracking down heterogeneity to design effective disease control strategies. As an example, let us consider a community where, prior to any intervention, 10% of the population has a 10-fold transmission risk and R 21 0 = 2 . 5 . Suppose that there exists some action (be it sustained chemoteraphy administration, improved hygiene or education; see e.g. Rollinson et al., 2013 ) by which a fraction η of the total population can be effectively removed from being susceptible to infection (independently of the risk group they belong to). Also, assume that the implementation of the interventions can be either evenly spread within the population or targeted (i.e. high-risk group members are prioritized). In the former case, the allocation of the population between the two risk groups does not change; in the latter, the high-risk share of the population is reduced to max (0 , ( f − η) / (1 − η)) . The top panels of Fig. 6 show that bringing R 21 0 below one, thus breaking transmission, would require treatment for 7% or 60% of the overall human population with targeted or untargeted interventions, respectively.
The dual case of heterogeneity associated with nonhomogeneous water contact patterns and/or snail distributions (source heterogeneity) may be seen as relatively less critical, in general, because the dilution effect induced by the availability of multiple water sources typically hinders long-term parasite establishment. However, one could still investigate possible differences between targeted (priority is given to the most used and/or snail-rich water sources) and homogeneous snail removal (be it by means of focal molluscicides or biological control; see King and Bertsch, 2015;Sokolow et al., 2015 ). We refer to the prototypical case of a single well-mixed human community with access to two water sources. As an example, let us consider a community in which, prior to any intervention, the overall snail population is evenly split between the two water points but access to water sources is heterogeneous, with 75% of contacts being made at one site (say site 1), and R 12 0 = 2 . 5 . In this case, interventions are assumed to be able to reduce snail abundance by a fraction ξ of the total population. The distribution of the snail population between the two water sources does not change in case of a homogeneous intervention. With a targeted intervention, instead, the fraction of snail present at the most used water source is reduced to max (0 , (n 1 − ξ ) / (1 − ξ )) . The bottom panels of Fig. 6 show that endemic transmission could be interrupted by removing about 34% of the overall snail population with a targeted intervention, instead of 60% in the case of a blanket intervention. Unsurprisingly, then, interventions to reduce schistosomiasis transmission are increasingly more efficient (and possibly cost-effective) if targeted to the most-at-risk groups of the human population (transmission hubs) or, to a lesser extent, to the most critical water sources (transmission hotspots; see Woolhouse et al., 1991Woolhouse et al., , 1998. Remarkably, both examples show that targeted actions can interrupt endemic transmission even if transmission hubs or hotspots cannot be treated completely. Therefore, not only does our study affirm that accounting for heterogeneity may be of paramount importance for a correct characterization of epidemicity thresholds, infection intensity patterns, and other quantities of chief epidemiological interest; it does also suggest that control measures designed under the assumption of homogeneity -of hosts' population abundances, transmission rates and/or mobility patterns -could be misleading, and that evaluating control strategies (and their implementation) using homogeneous modeling tools could severely over-or under-estimate their effectiveness. In real-world communities, many (if not all) of the sources of heterogeneity explored in this work are likely to be in play -possibly simultaneously. These considerations raise a crucial question about how the different types of heterogeneity can be quantified from field data. To answer this question, we note that proxies of transmission risk heterogeneity can be estimated via census data, which provide details on the demographic, social and economic structure of a community, as well as on its access to safe water and sanitation (see e.g. Mari et al., 2012Mari et al., , 2017. This type of information could be usefully complemented with in situ surveys to unveil (possibly season-dependent) water-related habits and behaviors (see e.g. Woolhouse et al., 1991Woolhouse et al., , 1998. Water-contact patterns can be assessed via interviews on the use of water sources, which in turn require a detailed mapping of water access points. The spatiotemporal distribution of the snail population in different water bodies can be evaluated via malacogical surveys and habitat suitability studies based on remote sensing and geospatial analysis (e.g. Brooker, 2007;Simoonga et al., 2009;Stensgaard et al., 2013;Pedersen et al., 2014;Perez-Saez et al., 2015 ). Human mobility patterns (daily commutes and seasonal migrations) can be inferred from ad hoc questionnaires or, on a larger scale, from the analysis of mobile phone data (see e.g. Wesolowski et al., 2012;Bengtsson et al., 2015;Finger et al., 2016;Ciddio et al., 2017;Mari et al., 2017 , for epidemiological applications). Tapping one or more of these sources of information could help identify transmission hubs and hotspots, thus paving the way for more effective intervention tactics. On a cautionary side, however, we note that targeting heterogeneity may but just one (yet possibly important) piece of the jigsaw puzzle that is schistosomiasis control, and that many other crucial factors (e.g. type(s) of intervention, timing, cost-effectiveness) must be properly taken into consideration as well (see e.g. Rollinson et al., 2013 ).
The modeling approach analyzed in this work is clearly not exempt from limitations. One additional source of heterogeneity that is neglected in Macdonald-like models for schistosomiasis transmission is the assumption that the average parasite burden in the human community (or in a relevant subset) is an informative measure of parasite distribution in the population. Individual variations in parasite load, which are often observed in field data, could be addressed by introducing a proper stratification of infection intensity ( Gurarie and King, 2014;Gurarie et al., 2010Gurarie et al., ,2016, at the expense of a relatively more complicated modeling frame-work. Also, human demographic dynamics, that are typically neglected in Macdonald-like models, can bear important implications for schistosomiasis dynamics and control. Human communities in developing countries may be far from demographic equilibrium ( Bongaarts, 2009 ), thus a static stratification of the population into age-groups may be not sufficient to fully represent the actual demographic dynamics. Given the much shorter timescales involved, spatiotemporal variations in snail demography (e.g. Perez-Saez et al., 2016;Gurarie et al., 2017 ) may perhaps play an even more important role in determining transmission patterns, especially in highly seasonal environments, where fluctuations of different variables (such as snail densities and human exposure rates) may be remarkably lagged. Indeed, when coupled to a realistic description of the snail intermediate hosts' population ecology, seasonal forcing can not only influence the likelihood of parasite establishment, but also induce a wide range of endemic dynamics, including periodic, quasi-periodic and chaotic transmission patterns ( Ciddio et al., 2015 ). Spatial coupling mechanisms are also underrepresented in this modeling framework, which, in its present form, can accommodate neither hydrological transport of the larval forms of the parasite ( Maszle et al., 1998;Lowe et al., 2005 ) nor snail dispersal ( Clennon et al., 2007;Wu et al., 2008 ). These shortcomings could be overcome at the expense of a more complex modeling approach describing several layers of spatial connectivity ( Gurarie and Seto, 2009;Perez-Saez et al., 2015;Ciddio et al., 2017 ), possibly also accounting for the geographical signatures of environmental fluctuations ( Mari et al., 2014 ). We finally note that a careful consideration of the spatial scale of analysis that is relevant for the application under study is crucial to determine what sources of heterogeneity are most likely to influence disease transmission.
Complemented with one or more of these possible extensions, the simple theoretical framework provided by Macdonaldlike models can indeed be turned into a powerful operative tool to gain insight into the role of heterogeneity in schistosomiasis transmission dynamics, as well to analyze real-world intervention strategies.