Schistosoma Transmission in a Dynamic Seasonal Environment and its Impact on the Effectiveness of Disease Control

Abstract Background A seasonal transmission environment including seasonal variation of snail population density and human-snail contact patterns can affect the dynamics of Schistosoma infection and the success of control interventions. In projecting control outcomes, conventional modeling approaches have often ignored seasonality by using simplified intermediate-host modeling, or by restricting seasonal effects through use of yearly averaging. Methods We used mathematical analysis and numerical simulation to estimate the impact of seasonality on disease dynamics and control outcomes, and to evaluate whether seasonal averaging or intermediate-host reduction can provide reliable predictions of control outcomes. We also examined whether seasonality could be used as leverage in creation of effective control strategies. Results We found models that used seasonal averaging could grossly overestimate infection burden and underestimate control outcomes in highly seasonal environments. We showed that proper intraseasonal timing of control measures could make marked improvement on the long-term burden reduction for Schistosoma transmission control, and we identified the optimal timing for each intervention. Seasonal snail control, implemented alone, was less effective than mass drug administration, but could provide additive impact in reaching control and elimination targets. Conclusions Seasonal variation makes Schistosoma transmission less sustainable and easier to control than predicted by earlier modeling studies.

Schistosomiasis, caused by Schistosoma species parasites, is highly prevalent in the tropics. Like many infectious diseases, Schistosoma transmission is highly seasonal with a timedependent transmission rate that varies during a year. Because schistosome parasites split their life cycle between humans (who harbor the adult forms of the parasite) and intermediate-host freshwater snails (which harbor and amplify the infectious larval forms), human water contact, area contamination with sewage, and the presence of competent snail species all play significant roles in defining local transmission [1]. The observed seasonal nature of Schistosoma transmission is often related to the variable abundance of intermediate host snail populations, which is linked to weather patterns (rainfall, temperature), and/or the resulting availability and frequency of human-snail contacts [2][3][4][5]. The mathematical model of Schistosoma transmission was first proposed by G. Macdonald [6] and later developed in different directions [7][8][9][10]. However, most mathematical models of Schistosoma transmission do not account for the effects of seasonality and this could significantly affect the accuracy of their predictions. A deeper understanding of the effects of transmission seasonality should provide additional insight into the optimal timing of interventions used in public health practice for schistosomiasis control. The current options for control include periodic mass administration of the drug praziquantel and, in some locations, control of intermediate host snails with the use of chemical molluscicides or environmental modification [11][12][13][14].
In mathematical models of Schistosoma transmission, one way to deal with seasonal variability has been to use seasonal averaging of environmental and behavioral inputs, making the resulting stationary Macdonald system amenable to direct analysis [15]. However, little effort has been put into assessing the effect of such seasonal averaging on the resulting predictions of disease-risk or control outcomes (some recent works have made such attempts, eg, [16]). In general, the use of averaging methods in nonstationary dynamical systems can be justified when such variability is marginal relative to the baseline mean state, or when the quantities of interest (eg, dynamic variables or outputs) depend linearly on variable inputs. Under these conditions, the averaged stationary model, with properly adjusted coefficients, is able to reproduce the approximate "mean" behavior of the nonstationary (time-periodic) system. However, for nonlinear disease transmission models in variable seasonal environments, the output histories can depart significantly from the expected "mean" behavior [16], because of highly influential seasonal effects on local transmission.
Another commonly used procedure in modeling vector-host transmission is a reduction of the coupled human-snail transmission cycle to a single host equation via quasiequilibration of the snail equation. This can be justified by the relatively short duration of snail lifespan, relative to human and worm, but requires more careful analysis, both for the basic Macdonald system and for its extensions.
In the present study, we aimed to perform a systematic assessment of seasonal variability and its implications for Schistosoma control predictions, including those for programs relying on mass drug administration (MDA) and/or molluscicide-based snail control to achieve their public health objectives.

METHODS
We modeled seasonal variability as seasonal changes of snail abundance, knowing that multiple environmental factors (precipitation, temperature, etc. [3][4][5][17][18][19][20][21]) can affect snail population density. In our modelling analysis, we used 2 types of snail population dynamics: (1) a prescribed periodic snail density function and (2) a dynamic snail population model based on the underlying periodic carrying capacity function, combined with estimates of snail reproduction rate and mortality.
The key dimensionless input parameters of the model were the basic reproduction number (R 0 ), and the amplitude of seasonal variability. We used type 1 snail population models to explore the parameter space of the system and identify regions of sustained infection versus elimination. For control analysis, however, we employed the dynamic snail model (type 2) to account for abrupt changes in population due to molluscicide. We then explored different periodicity patterns to approximate seasonal snail population dynamics, as informed by empirical data from past field studies on snail abundance [4,5,17,18,22,23]. Details of the models are provided in Supplementary Material.

Modeling of a Stationary Transmission Environment and Derivation of R 0
The basic Macdonald model for Schistosoma transmission combines 2 variables: mean worm burden (MWB) w (t) of human hosts, and infected snail prevalence y (t), coupled by a system of differential equations: Here, transmission coefficient A represents mean rate of worm accumulation per human host and B the force of snail infection [15,24,25]. Coefficient A is proportional to the snail density coefficient, N, coefficient Bdepends on the human population size H, and its infectivity (egg shedding by mated worm pairs in the human population). Both coefficients are proportional to human-sail contact rates. Coefficients γ and ν are the natural mortality rates for adult worms and snails. We took γ = 0.25/ year and allowed ν to vary from 6 to 2 per year (equivalent to a 2-month to half-year lifespan) [26,27].
The system (Equation 1) can be rescaled in dimensionless form, given by a single parameter R 0 = A B γ ν (see Supplementary Material). Parameter R 0 measures the intensity of transmission in a hypothetical stationary Macdonald system. Namely, it has stable endemic state (w * , y * ) > 0 when R 0 > 1, and a stable infection-free equilibrium (w * = y * = 0) when R 0 < 1 . The analysis in the current paper mostly employed the basic Macdonald system; however, several modifications are discussed in the Supplementary Material.
The nonstationary (seasonal) transmission environment, the subject of our present work, brings another important parameter-amplitude of seasonal variability [28]. Such seasonality can arise from multiple sources, that is varying snail populations or seasonally varying human behavior (water exposure/ environmental contamination). Here, we focused on seasonal variation in snail population dynamics and its implications for sustained transmission and control. The potential influence of other variable inputs is discussed in the Supplementary Material.

Incorporating a Seasonal Snail Environment
Multiple environmental factors can affect snail population dynamics (reproduction, mortality/survival, development, and infection). They include temperature, rainfall ( Figure 1A), food resources, agriculture, snail population predation, and environmental chemistry [3,4,19,[29][30][31][32][33]. We did not include all such factors explicitly, but considered their effects on snail dynamic in 2 ways: (1) a prescribed periodic snail population function N(t, a) of seasonal amplitude, a; and (2) a logistic snail population model driven by an underlying carrying capacity function K(t, a).
In both cases, the variability of N and K was determined by their amplitude parameter 0 ≤ a < 1, where a = 0 corresponded to a stationary population value, while a larger a > 0 marked the departure from this value. For our analysis, all population density functions were dimensionless, having been rescaled relative to a putative mean population density N * .
The simplest mathematical form of periodic functions (N, K) is trigonometric, referred to in this paper as type-I seasonality. The distinguishing features of type I are evenly distributed high and low  [3,22]. B, Examples of type-I (left) and type-II peak (right) seasonality with amplitude parameter a. The type-I seasonality was modeled by trigonometric function, 1 + a cos (2 π t), and type-II seasonality was modeled by an elliptic theta function of amplitude 0 ≤ a < 1, 1 + 2 ∞ n=1 a n 2 cos (2π n t). At small amplitude, a, both types are approximately equal, because higher-order Fourier modes become negligible. But as amplitude increases, they depart significantly in their variability (finite for type I and unlimited for type II). C, Seasonal average of mean worm burden (MWB), w (a) = w * (t, a) as a function of amplitude, a, for human-snail Macdonald systems Equation (4) for 3 values of basic reproduction number R 0 (transmission intensity), R 0 = 1.5; 2; 3. Left, results for a type-I trigonometric N (t, a) model; right, results for a type-II peak N(t, a) model. The curves indicate that for a lower R 0 and a higher seasonal amplitude, transmission becomes unsustainable in both type-I and type-II models. Not shown, these curves also depend on snail mortality, which in the case shown was ν = 4. seasons, where the snail population varies about its mean value (= 1) in equal proportions and evenly spaced across high and low seasons. This type-I system has limited range of seasonal variability (amplitude or variance).
However, in many cases, more extreme patterns of seasonality arise, for example, a short rainy season with high snail density, followed by longer dry periods of low density, or multiple seasonal peaks [22]. A convenient functional form for such variability can be described via an elliptic theta-function of amplitude (0 ≤ a < 1), Such N becomes highly concentrated as a approaches 1 (Dirac delta-function) ( Figure 1B). We refer to it as type-II (or peak) seasonality. The (nonlinear) amplitude of type-II (Equation 3) has different meaning than type I. Both types are approximately equal at small a, but as amplitude increases, they depart significantly in their variability. Data on snail abundance supports both patterns [4,5,17,18,22,23]. Of course, they represent a crude approximation of such variability. Their simple analytic form is convenient for analysis and numeric simulation.

The Macdonald Model System With Seasonal Snail Variation
The rescaled seasonal Macdonald system with a prescribed, variable snail population has the form These equations are derived similar to the stationary case, but with modifications to account for nonstationary snail density, N (t, a). Specifically, the infected snail prevalence variable 0 < y (t) < 1 of the stationary model is replaced by an infected snail density (0 < y (t) < N (t, a)), and the constant snail mortality ν is replaced by time-dependent function, µ (t).
Here β is maximal reproduction rate, K (t, a) is seasonal carrying capacity, and ν is natural snail mortality [22]. Snail dynamic (Equation 5) will automatically reproduce time-dependent mortality, µ(t).
The snail population in Equation 5 is thus decoupled from the standard Macdonald system of Equation 4. We note that snail infection typically has only a marginal effect on its overall population rates for reproduction/growth/mortality due to relatively low levels of patent (shedding) snails in natural environments [34,35].
As with N (t, a), the periodic function K (t, a) can be either a trigonometric type (type I) or a peak type (type II). The logistic Equation 5 with periodic K (t, a) has a stable periodic solution N * (t, a), which plays the role of the prescribed function N (t, a) in the Macdonald case (Equation 4).

Macdonald System
A single MDA session at time T in a host community with MWB w (t)was simulated in our model as an instantaneous where the treatment effectiveness constant, ε H , combines antiparasite drug efficacy (the fraction of killed worms), ε, and the population coverage fraction 0 < f < 1, that is ε H = εf . Such formulation implies each MDA session draws a random host pool (fraction f) for treatment. Another possibility is systematic noncompliance whereby the same population is excluded from treatment. The latter requires an extension of the Macdonald system discussed in Supplementary Material A16.
Instantaneous reduction of worm burden post MDA is due to short half-life (hours) of praziquantel relative to time scales of infection history (months to years). It is implemented via reinitialization of variable w (t) at time T. Snail control via molluscicide (at time T) was implemented in a similar fashion as an instantaneous event, whereby dynamic variables {N (t, a) , y (t)} are dropped by factor 0 < ε S < 1 (reflecting the snail killing fraction), which depends on the efficacy of molluscicide application,

Seasonal Snail Variability and its Effects on Endemicity of Human Infection
The nonstationary Macdonald system with seasonally varying snail population, N (t, a), has 2 dimensionless parameters: R 0 (intensity of transmission) and amplitude (0 ≤ a < 1). We want to study their effect on sustainability of transmission. The role of stationary endemic equilibria is demonstrated here by timeperiodic solutions{w * (t, R 0 , a) , y * (t, R 0 , a)}. Their stability types and dependence on parameters (R 0 , a) space is studied via seasonal average of the MWB-variable w (t, R 0 , a) →w (a, R 0 ). The stationary case corresponds to (a = 0). All results below are obtained via numeric simulations of the time-periodic dynamical Macdonald system. Key findings are listed below (see Supplementary Material for details): 1. w (a, R 0 ) is a decreasing function of a, with maximum w max = 1 − 1/R 0 , attained at a = 0 (ie, stationary case), and minimal value at large values (a ≈ 1) ( Figure 1C). The drop of w (a, R 0 ) along the a-axis depends on the transmission intensity (R 0 ). It drops faster for smaller R 0 , and can drop to 0, which means infection becomes unsustainable. Figure 1C shows functions w (a, R 0 ) for both types of seasonality and 3 values of R 0 . For type I we observed a 15% drop of (w (1, R 0 ) /w (0, R 0 )) for R 0 = 3, and a larger 35% drop and 75% drop for R 0 = 2 and R 0 = 1.5, respectively. This R 0 /amplitude effect was even more pronounced for type-II seasonality, where all periodic stable solutions become unsustainable at sufficiently large a.

The periodic Macdonald system (Equation 4) behaves qualitatively similar to a stationary Macdonald case (Equation 1).
Specifically, it has infection-free equilibrium (w = 0) and periodic endemic states (w (t, R 0 , a)), whose stability types are determined by parameters (R 0 , a). Theoretically, the (R 0 ,a) space can be divided into 2 regions of possible long-term outcomes: a stable infection-free sector and a stable endemic state sector. But in practice, it is difficult to find the boundary numerically. In Figure 2 we show 2 types of parameter space analysis based on infection-free equilibrium and periodic endemic state, respectively. In both cases, increased seasonal amplitude makes infection less sustainable; transitioning from stable endemic to stable infection-free is gradual and difficult to assess computationally (see Supplementary Material for details).
In brief, this suggests that seasonal variability of snail numbers could make local transmission and persistence of infection potentially less sustainable and, thus, ultimately easier to control.

Optimal MDA Timing
We next explored the impacts of 2 different types of control interventions, MDA and molluscicide. The former has no direct effect on snail populations, while the latter targets snails specifically. This analysis required use of the dynamic snail population model, that is the coupled system of Equations 4 and 5. The key input here was the seasonal carrying capacity function, K (t, a), represented by the type-I or type-II periodic functions.
In all cases, we used the endemic periodic state {w * (t, a) , y * (t, a)} to initialize the system, and ran a 6-year control program, as suggested by the World Health Organization (WHO) [14,36]. Specifically, annual MDA has an effective parameter ε H = ε f , which combines drug efficacy (ε = 0.70 − 0.85) and the population coverage fraction ( f ). Outcomes of such intervention depend foremost on transmission intensity, R 0 . Here we took an intermediate value (R 0 = 3) [36], and strong seasonality, amplitude a = 0.9 for type I, and a = 0.6for type II. Other examples are discussed in Supplementary Figures 12-17).
We found that different choices for intraseasonal MDA timing during the annual cycle (time =0 ≤ τ < 1) can affect the MWB reduction. The upper panels of Figure 3 demonstrate this effect by comparing MDA histories implemented at  ). Here we used the MDA efficacy ε H = 50%, corresponding to drug efficacy (70%-85% reduction in MWB) combined with the MDA treatment population coverage fraction (60%-70%). The upper panels show 6-year histories with MDA given at the start of each season τ = 0 (blue), or at midseason τ = 1/2 (yellow). The lower panels show the 5-year seasonal averages for MWB among local humans when different seasonal timings (0 ≤ τ < 1) were used for implementation. For each seasonal timing, τ = 0, 1 10 , 1 4 , 1 2 , 3 4 , 9 10 , we averaged MWB over the proper 1-year time interval, ie, over the period (τ , τ + 1), for successive 6 years. Qualitatively, these results look similar to R 0 = 6 (high transmission intensity) case discussed in Supplementary Figure 14, but the seasonal difference is less significant. τ = 0 (the season's start, set up at the peak of snail population) versus τ = 0.5 (midseason). To estimate the optimal MDA timing, that is the greatest MWB reduction by year 6, we ran multiple simulations with different choices (0 < τ < 1). In each case, we averaged the resulting MWB-solution {w (t)} over the intertreatment time interval, [τ , τ + 1].
In the lower panels of Figure 3, we show the resulting average MWB values over a 6-year control history for different choices of MDA implementation time, τ . The difference between the "best" and "worst" timing, was marginal-a 2.5% MWB reduction by year 6. Such an effect, however, could be more significant in areas with high transmission intensity, for example, R 0 = 6 (Supplementary Table 4). The optimal timing for MDA administration varied in the range, 1/2 < τ < 3/4, that is close to the midseason for high amplitude (a) but shifted toward τ ≈ 3/4 for moderate a. The pattern was similar for both types of seasonality. The results are summarized in Table 1, which shows, for each τ tested, the post intervention average MWB (w (Y 6 , τ ) /w (Y 0 , τ )) endemic values.

Optimal Timing for Molluscicide Application and for Integrated Control Strategies
We next used our dynamic snail models to estimate the optimal timing for molluscicide-based snail control, and for an integrated strategy of MDA plus snail control.
We used the same R 0 and a values as above. Snail mortality and maximal growth were fixed at ν = 4/year, β = 20/year, broadly consistent with some estimated parameters [22]. Figure 4 shows a projected 6-year history for a molluscicide-only control program, comparing type-I and type-II seasonality simulations. Two choices of molluscicide timing are compared, τ = 0 (blue) and τ = 1/2 (yellow). This effect of molluscicide on human infection (8%-38% reduction in MWB; Table 2) was less significant than drug treatment (82%-98% reduction in MWB; Table 1). Nevertheless, if snail control was accompanied by introduction of optimally timed MDA, human MWB (Table 3) and infected snail density y (t) (not shown) have approached near-elimination state after a 6-year program.
The optimal timing of molluscicide implementation for maximal MWB reduction ( Table 2) fell near the start of the season (τ ≈ 0) where carrying capacity function, K (t, a), and snail population density, N (t, a), approach maximal values. The overall progress measured as above by w (Y 6 , τ ) /w (Y 0 , τ ) varied in the range 8%-38%, depending on the type and amplitude of seasonality and the efficiency, ε S , of molluscicide. We also observed that the effect of optimal timing was much more significant for molluscicide than for MDA (Table 2).
A combined MDA plus molluscicide strategy was next simulated for 3 choices of molluscicide efficacy, ε S = {50%, 70%, 90%}, and 3 choices for MDA efficacy, ε H = {50%, 55%, 70%}. For this analysis, each intervention was to be given at its own optimal timing, τ , as suggested by the results in Table 1 and Table 2, namely τ = 0.5 − 0.75 for MDA and τ = 0 for molluscicide. We then compared 2 optimized strategies, MDA alone versus MDA plus molluscicide in Table 3. This indicates that integrated control can bring a significant improvement above an MDA-only strategy when each intervention is properly timed. Depending on seasonal amplitude, Table 1 Data are percent reduction in MWB.

. Effect of Annual Mass Drug Administration (MDA) in Terms of Mean Worm Burden (MWB) Percent Reduction Among Local Human Populations
We simulated a 6-year MDA regimen for a Macdonald-type model system having dynamic snail populations with seasonal carrying capacity function K (t, a) of trigonometric types I or peak type II. Impact was measured in terms of relative MWB reduction for local residents (year 6 over endemic year 0) using seasonal average values w (Y6, τ ) /w (Y0, τ ) of MWB. We examined several different values of intraseasonal timing for MDA, τ = 0, 1 10 , 1 4 , 1 2 , 3 4 , 9 10 , as fractions of the seasonal cycles, 3 possible levels of MDA efficacy (instantaneous reduction of mean worm burden) εH = {50%, 55%, 70%} as the combined effect of drug efficacy (fraction of killed worms) and population coverage fraction, and 2 values of seasonal amplitude, a. The optimal timing (highlighted in bold) varied in the range 1 2 < τ < 3 4. The optimum was close to mid-season τ ≈ 1 2 in the presence of strong seasonal amplitude (a) and shifted toward τ ≈ 3 the type of seasonality, transmission intensity R 0 , molluscicide clearing efficiency ε S , and MDA efficacyε H , the year 6 worm burden could be brought to elimination, relative to the baseline (endemic) values.

DISCUSSION
Conventional modeling approaches to infection transmission dynamics in a seasonally varying environment often employ time-averaged stationary models [15]. Such approximation can be justified when seasonal variability is low relative to the mean state but may not hold in general. Here, we studied the effects of seasonal snail population using Macdonald-type model with 2 basic patterns of variability, trigonometric (type I), and peak (type II). Both have an explicit mathematical form, but our analysis has employed primarily numeric simulations of the appropriate dynamical systems.
The key inputs for our analysis are dimensionless parameters: (1) the basic reproduction number of the stationary (seasonal-mean) system, R 0 (intensity of transmission); and (2) the amplitude (a) of seasonal snail variability. We addressed several questions of the combined effect of (R 0 , a) on periodic and stationary patterns of human and snail infection levels, and their seasonal averages, exemplified by MWB function w (R 0 , a). As function of amplitude, a, it typically maintained a plateau region at small or moderate a values, so seasonal averaging for Macdonald system would still be approximately valid. However, it would fail for large amplitude, where w (R 0 , a) underwent a rapid decay at high amplitude values, in some cases to w = 0 (elimination). Along these lines, we identified specific parameter ranges where transmission becomes unsustainable. Overall, we found increased seasonality makes infection transmission and endemicity less sustainable, hence more amenable to control interventions, as compared to areas with unchanging (stationary) annual transmission. Future studies could explore whether such conclusions could hold for other environmentally and seasonally mediated infectious diseases,  Here we used molluscicide efficacy (percent of killed snails) ε S = 70%. The 2 colors correspond to different seasonal timing of molluscicide application: at the start of the season, τ = 0 (blue) or at midseason, τ = 1 2 (yellow). The effect of seasonal timing on long-term patterns of transmission (6-year history) can be significant: τ = 0 implementation ultimately gives higher worm-burden reduction among local humans than does τ = 1 2 implementation. Qualitatively, these results look similar to R 0 = 6 (high transmission intensity) case discussed in Supplementary  Figure 15, but the seasonal difference is less significant.
like malaria. Most of our analysis was done using the basic Macdonald system, but some modifications were extended to evaluate other transmission models, including a Macdonald system with worm mating [25] and the stratified worm burden system [37,38] with comparable results (see Supplementary Material for details).
Seasonality can affect snail population dynamics and human behavior (environmental exposure and contact rates). The 2 factors are sometimes combined into a single seasonal transmission rate [16,39]. There are, however, significant differences between them (Supplementary Material), and the current analysis has specifically focused on seasonal snail dynamics.
We have applied our models to leverage seasonality for the optimal timing of repeated control interventions (including annual MDA, molluscicide application, or integrated intervention strategies that combine the 2 approaches) to predict which approach could achieve maximal reductions in parasite burden over a limited time span (6-10 years). The optimal intraseasonal timing (0 < τ < 1) varied for different interventions; the season start was defined by maximal snail density. For MDA alone, we A 6-year control program was simulated using a Macdonald-type model system having a dynamic snail population seasonality with carrying capacity function K (t, a) of trigonometric type I or peak type II structure. As for Table 1, progress was measured by relative seasonal average reduction of MWB w (Y6, τ ) /w (Y0, τ ). We examined several possible intraseasonal timings for molluscicide application,τ = 0, 1 10 , 1 4 , 1 2 , 3 4 , 9 10 , as fractions of the seasonal cycle, 3 possible levels of molluscicide efficacy (percent of killed snails) ε S = {50%, 70%, 90%}, and 2 choices of the seasonal amplitude parameter (moderate or high). Optimal timing for molluscicide control in all cases was τ = 0 (highlighted in bold), ie, at the start/end of the season, when the snail population or its carrying capacity reached its maximum. Data are percent reduction in mean worm burden (MWB).
As in Tables 1 and 2, progress was measured by relative seasonal average reduction of MWB w (Y6, τ ) /w (Y0, τ ). A 6-year control program was simulated for a Macdonald-type model system having dynamic snail populations (K (t, a) of type-I or type-II), 3 levels of MDA efficacy εH = {50%, 55%, 70%}, 3 levels of molluscicide efficacy ε S = {50%, 70%, 90%}, and 2 different amplitudes of seasonality. The impact for MDA-only was estimated at its optimal time of delivery, as suggested by Table 1. The impact of the MDA + snail strategy was based on separate delivery at their individual optimal timings, τ , as suggested by Table 1 and Table 2, namely τ = 1 2 or 3 4 for MDA, and τ = 0 for molluscicide application.
found the optimal timing to be in the range, 1/2 < τ < 3/4, near the minimal snail density (or its carrying capacity). For molluscicide-based control, the optimal timing of delivery was closer to the start of season, that is at the snail population peak. The overall effect of mollusciciding on human worm burden reduction was less significant than MDA at short duration, even at high snail killing efficacy. But the effect of optimal timing (τ ) was more pronounced for molluscicide than for MDA. The combined strategy (MDA plus molluscicide), with proper timing for each, was found to provide an enhanced reduction of MWB compared to MDA alone. This suggests that addition of snail control can add substantially to the program impact, in particular at high transmission intensity where MDA alone may not be sufficient (Supplementary Tables 4-6).
In general, control outcomes depend on transmission intensity, strength and patterns of seasonality, drug and molluscicide efficacy, and control duration. A 6-year program is an intermediate-range based on WHO guidelines [14,36]. In some cases, it could lead to near-elimination after a 6-year period. In other cases, additional treatment or an integrated strategy combining MDA with molluscicide is needed. These results are consistent with the results from earlier modeling papers [36,40]. Another important factor is the simulation of the MDA delivery for a prescribed population fraction. For example, it could be a randomly drawn pool of residents or involve a systematically noncompliant pool. We observed that systematic noncompliance could significantly slow the progress of MDA program (see Supplementary Material for details). Identifying robust and optimal combination strategies for schistosomiasis control under these various conditions is paramount for achieving and maintaining WHO schistosomiasis elimination goals. Future work should include the development of systematic optimization procedures for identifying optimal control interventions across transmission settings.
We believe the basic conclusions on seasonal timing are robust and would hold in different environments, for different transmission intensity, R 0 , and for different schedules of MDA delivery. Furthermore, we expect they would hold for other transmission models, and with more realistic patterns of snail population dynamics driven by temperature and rainfall data [3,22,35,41].

Supplementary Data
Supplementary materials are available at The Journal of Infectious Diseases online. Consisting of data provided by the authors to benefit the reader, the posted materials are not copyedited and are the sole responsibility of the authors, so questions or comments should be addressed to the corresponding author.