Macro-parasite transmission in dynamic seasonal environment: Basic Reproductive Number, endemicity, and control

Seasonality of transmission environment, which includes snail populations and habitats, or human-snail contact patterns, can affect the dynamics of schistosomiasis infection, and control outcomes. Conventional modeling approaches often ignore or oversimplify it by applying ‘seasonal mean’ formulation. Mathematically, such ‘averaging’ is justified when model outputs/quantities of interest depend linearly on input variables. That is not generally the case for macroparasite transmission models, where model outputs are nonlinear functions of seasonality fashion. Another commonly used approach for Schistosomiasis modeling is a reduction of coupled human-snail system to a single ’human equation’, via quasi-stationary snail (intermediate host) dynamics. The basic questions arising from these approaches are whether such ‘seasonal averaging’ and ‘intermediate host reduction’ are suitable for highly variable/seasonal environments, and what implications these methods have on models’ predictive potential of control interventions. Here we address these questions by using a combination of mathematical analysis and numerical simulation of two commonly used models for macroparasite transmission, MacDonald (MWB), and stratified worm burden (SWB) snail-human systems. We showed that predictions from ‘seasonal averaging’ models can depart significantly from those of quasi-stationary models. Typically, seasonality would lower endemicity and sustained infection, vs. stationary system with comparable transmission inputs. Furthermore, discrepancies between the two models (‘seasonal’ and its ‘stationary mean’) increase with amplitude (or variance) of seasonality. So sufficiently high variability can render infection unsustainable. Similar discrepancies were observed between coupled and reduced ’single host’ models, with reduced model overpredicting sustained endemicity. Seasonal variability of transmission raises the question of optimal control timing. Using dynamic simulation, we show that optimal timing of repeated MDA is about half season past the snail peak, where snail population attains its minimal value. Compared to sub-optimal timing, such strategy can reduce human worm burden by factor 2 after 5-6 rounds of MDA. We also extended our models for dynamic snail populations, which allowed us to study the effect of repeated molluscicide, or combined strategy (MDA + molluscicide). The optimal time for molluscicide alone is the end or the start of season, and combined strategy can give additional reduction, and in some cases lead to elimination.


Introduction
Seasonal transmission of an infectious disease is a common phenomenon, defined by time dependent transmission rates that vary over the season. Schistosomiasis is a case example where the seasonality is related to snail populations linked to weather patterns (rainfall, temperature), and/or behavioral factors (human-snail contacts). But mathematical models of Schistosome transmission rarely account for seasonality, with implication to accuracy of model prediction and opportunity to gain insight into public health practice, like optimal timing of interventions.
Basic reproduction number ( 0 R ) is a fundamental concept widely used in population biology and infectious disease modeling. In the context of infectious diseases, 0 R is used to assess the magnitude of disease outbreak, spread or endemic (equilibrium) level [ 1 ]. Its applications range from modeling of human-to-human communicable diseases to vector-mediated systems for macroparasite diseases, such as Ross-MacDonald models of schistosomiasis (see e.g. 2 , 3 , 4 , 5 , 6 ). Though 0 R was first explicitly derived for simple single population models, this concept was further extended to more complex metapopulation systems, e.g. spatially connected environment with multiple host sites (see e.g. 7 , 8 , 9 , 10,11 ).
Mathematically, 0 R is most appropriate for stationary environment where disease transmission parameters are kept fixed. However, in many cases, such environment is highly variable, e.g. seasonal change of vector population or contact patterns, as well as abrupt changes caused by control interventions.
One way to incorporate such variability in mathematical models is through appropriate 'seasonal averaging' of environment and/or behavioral inputs. However, little effort was spent to assess the effect of 'seasonal averaging' for analysis of disease predictions, and control outcomes. In general, the use of averaging methods in non-stationary dynamical systems can be justified when variability is marginal relative to baseline mean state or when quantities of interest (e.g. dynamic variables or outputs) depend linearly on variable inputs. Under such conditions, the 'averaged' stationary model with properly adjusted coefficients is able to reproduce approximate 'mean' behavior of the non-stationary (time periodic) system. However, neither of these conditions holds for nonlinear disease transmission models in variable seasonal environment. The outputs of such system are nonlinear functions of its input variables, and they can depart significantly from the expected 'mean' behavior.
Besides seasonal averaging, another commonly used approximation in modeling vector-host transmission is a reduction of the coupled system to a single host equation, via quasi-equilibration of fast vector variables. Such reduction is commonly justified by distinct time scales: slow 'host-parasite dynamics' vs. fast 'vector-intermediate host '. Indeed, in the context of schistosomiasis time scales for snail dynamics vary from weeks-to-months, while human and worm scales are order of magnitude slower (years).
Reduced models are convenient for theoretical analysis and for numeric simulations (see e.g. 2 , 3 ). Such reduction procedures however, require more careful analysis beyond simple 'time-scale' heuristics, to assess possible departure of the reduced or averaged models ('single-host' or 'seasonal-mean') from their fully coupled, non-stationary precursor.
Here we assess the effect of seasonal averaging and intermediate host reduction for schistosomiasis transmission in seasonal environment, and their implications on predicting control outcomes (MDA). On human side, we shall mostly use mean worm burden (MWB) MacDonald formulation ( 4 , 5-7 ), but make a few comment on other possible approaches, like stratified worm burden system ( [8][9][10]. Seasonal transmission environment can appear as variable snail population and/or variable human-snail contacts. In this paper, we mainly focus on seasonal snail variability (their population density) that determines their infection level, and snail-to-human transmission rate. Mathematically, different functions can be used to account for variable snail population, under different environmental conditions(precipitation, temperature, et al) (see e.g. 11,12 , 13 , 14 , 15 , 16 ). The two commonly observed patterns include: (i) moderate snail variability about 'mean value' with approximately equal high and low seasons, (ii) long dry season (with low / zero snail density and transmission) followed by short rainy season, with bursting snail population and intensified transmission rate.
We shall employ two mathematical functions for such seasonal patterns, namely (i) trigonometric function ( ) 1 cos 2 a t π + with amplitude 0 1 a < < , for alternating high/low seasons of equal duration, (ii) concentrated peak density over low-level background with prescribed seasonal mean value =1. The latter also have a version of amplitude parameter ( 0 1 a < < ), but its meaning and implication are different.
Earlier theoretical papers ( 17 , 18 ) studied seasonally varying (periodic) human-snail contact patterns and their impact on schistosomiasis endemic state and persistence. In particular, 18 estimated endemicity thresholds for periodically varying MacDonald system. Here we change the focus from 'seasonal contact' to 'seasonal snail population', and extend the scope of their analysis by exploring the implications of seasonality for sustained infection, and control interventions, MDA and molluscicide.
The bulk of our analysis employs prescribed snail population (density) function ( ) N t . The study of molluscicide control however, requires a dynamic snail population that could respond to an abrupt change, environmental and/or human-made interventions. Such resource-driven models of snail population biology were developed in several papers (e.g. 15 , 16 , 19 ). These works however, did not apply dynamic snail models for transmission/control analysis. Our paper attempts to fill this gap. Alongside prescribed function ( ) N t we develop model of resource driven seasonal snail, and apply it to study molluscicide.
Two basic parameters employed in our analysis are basic reproduction number 0 R (of 'mean' stationary system), and amplitude (a) of seasonal variability.
By combining mathematical analysis with numeric simulations, we address several questions: Parameter 0 R is responsible for endemic equilibria of (1.2) and for endemic to zero transition. A simple illustration is herd immunity, vaccine coverage fraction ( ) to prevent an outbreak.
Second order system (1.2) can be reduced by replacing 'fast' snail equation is with its quasi-equilibrium value (see Box 1) due to' short' snail time scales (months) compared to 'slow' human-worm dynamics (years). The reduced MacDonald system is a single equation Both full and reduced systems (1.2)-(1.3) share endemic equilibrium in a stationary environment.
Significant differences however, could arise is seasonal environment.
MacDonald system with mating and aggregation. Many conventional approaches, going back to the original MacDonald paper 4 , assume worm burden is distributed in host population according to negative binomial (NB) pattern with mean value , and aggregation parameter 0 k > ( 5,7,20,21 In stationary environment such (H-S) system has stable endemic equilibrium, Compared to (1.2) modified system (1.5) depends on 2 dimensionless parameters, transmission , rather than single 0 R . There are other fundamental differences between simple MacDonald (1.2), and its modified version (1.5). The former has two equilibria (zero-endemic), whose stability types are determined by 0 above critical value c R . Figure

Seasonal environment
Non-stationary transmission dynamics can arise from multiple sources, including natural (seasonal snail populations), behavioral (human-snail exposure/contamination), or repeated interventions, like MDA or molluscicide. In this section we focus on variable snail population, given by function N(t) (see e.g. 19 , 15 ). There are different types of snail environment, that make time dependent snail population ( ) N t . We do not attempt to classify and study all such patterns of natural variability, but confine our analysis to 2 typical periodic functions I.
. CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. not certified by peer review)

II.
Peak-type ( ) , N t a of amplitude 0 1 a ≤ < (Figure 1), defined via elliptic theta function. Amplitude a has different meaning now. Large a (close to 1) features sharp seasonal peak of short duration about t=0, followed by long dry period ( ( ) 0 N t ≈ ) over the remaining 'dry' season.
Data from different studies on snail abundance support such patterns, e.g. ( 13 , 19 22 ). In both cases, ( ) , N t a is rescaled so that its seasonal mean, ( ) . CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. not certified by peer review) (which was The copyright holder for this preprint this version posted February 3, 2020.  Figure 1). Some modifications are needed to accommodate variable snail population, namely and susceptible pool is given by clipped function (ii) fixed snail mortality (= rescaled rate 1 of (1.2)) is replaced by time dependent rate function . CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. not certified by peer review) (which was The copyright holder for this preprint this version posted February 3, 2020. . https://doi.org/10.1101/19012245 doi: medRxiv preprint Its reduced version is given by The key inputs for (1.7)-(1.8) is seasonal population density ( ) N t , and stationary-mean 0 R of rescaled system (Box 1). The role of conventional endemic equilibria ( sometimes combined into a single product-type term (see e.g. 17 , 18 ). While superficially 'seasonal snail' and 'seasonal contact' may look similar, there are significant differences in their equations, e.g. excess snail mortality in (1.7), (1.8); their dynamic responses may differ. In this paper we shall confine our analysis to seasonal snail population, and seasonal effect of periodic repeated MDA interventions.

Dynamic snail population
For most parts of our analysis, we use prescribed seasonal snail function ( ) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. not certified by peer review) (which was The copyright holder for this preprint this version posted February 3, 2020. . https://doi.org/10.1101/19012245 doi: medRxiv preprint growth model with prescribed maximal reproduction rate β , seasonal carrying capacity (CC) ( ) and mortality ν , Periodic CC-function ( ) , K t a is taken as trigonometric (type I), peak (type II). Equation (1.9) has stable periodic solution ( ) , N t a , which plays the role of the above prescribed type I-II N. We further assume snail infection has marginal effect on its growth/ mortality (due to small patency conversion fraction of snails). So equation (

Seasonal variability and sustained infection
We list a few key results of our analysis (further details can be found in Supplement).

i)
Stationary full (host-vector) and reduced (single host) MacDonald systems share the same equilibria and 0 R , but when perturbed from equilibrium (e.g. via MDA) their relaxation patterns . CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. not certified by peer review) (which was The copyright holder for this preprint this version posted February 3, 2020. . https://doi.org/10.1101/19012245 doi: medRxiv preprint to equilibrium differ, with reduced system relaxing faster (Supplement Figure A). Higher relaxation rate (rebound) would imply less efficient MDA control in reduced system compared to full MacDonald host-vector system. ii) The role of stationary equilibria in variable environment is played by periodic (dynamic) , , , w t a y t a , see Figure 2 (more details shown in Supplement Figure   Overall, we find seasonal snail variability makes infection less sustainable (c.f. 18 ), and thus easier to control. The reduced host only system overestimates sustainability and control outcomes for MacDonald system. . CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. not certified by peer review) (which was The copyright holder for this preprint this version posted February 3, 2020. . https://doi.org/10.1101/19012245 doi: medRxiv preprint

MDA control in stationary environments
Periodic MDA in stationary environment for simple MacDonald system An MDA session at time T is represented in our setup as instantaneous reduction of mean , where constant f ε combines drug efficacy ε (fraction of surviving worms), and population coverage fraction 0 The latter applies to population groups (e.g. worm strata), or the entire community (MacDonald MWB).
One can think of such MDA-event as a sharp spike of worm mortality over short duration. Formally represented by Dirac delta function.
Each MDA-event is followed by rebound -relaxation towards endemic state (equilibrium, periodic, et al).
Such relaxation patterns can differ in complete (2D) vs. reduced MacDonald (Supplement Figure A). A regularly spaced sequence of MDA-events creates yet another type of periodic variability in dynamic system. This time periodicity is affected by worm mortality, rather than snail population growth/decay. So ( ) t γ becomes a periodic function with sharp (Dirac delta) spikes, or their finite approximation. As above, one can ask whether such periodic mortality ( ) t γ can be approximated by its 'mean' value, ( ) Here the natural 0 R (Box 1) is replaced by effective reduced value due to (1.10) ( ) ( ) In particular, we ask whether reducing 1 ef R < (via suitable combination of MDA frequency and efficacy-coverage) could lead to elimination. The condition for elimination in the stationary 'mean' system is CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. not certified by peer review) (which was The copyright holder for this preprint this version posted February 3, 2020. To test the validity of MDA averaging, we run numeric simulations of two model, periodic MDA, and its 'average model" with enhanced worm mortality (1.10). The results (Supplement Figure E) show elimination path for the mean-MDA system, predicted by 1 ef R < , while the exact MDA-response curve is locked in a periodic (limit) cycle pattern, even above critical frequency ( c T T < ). Such limit-cycle response patterns arise in many transmission models, including SWB (see 23 ), they can frustrate target reduction goals via repeated MDA regimen (see 24 ).
Parameter space analysis for MacDonald system with mating.  Figure F). For detailed explanation, see ( 1 , 9 ). Another important feature of mated MacDonald system, shared by SWB 9 , is their dependence on two dimensionless parameters, e.g.
The bistable nature (breakpoint) of mated MacDonald system has implications for MDA control.  Figure G). In some cases (lower 'snail/human' ratio), it goes to elimination after finite number of MDA cycles (though total duration could vary). In other cases (high 'human/snail' ratio), the system is locked in a limit cycle, due to post-MDA rebound.
. CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. not certified by peer review) (which was The copyright holder for this preprint this version posted February 3, 2020.  .5 τ = (middle). It raises the question of optimal choice of MDA timing, to achieve maximal burden reduction by the end of the program (Y6).
. CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. not certified by peer review) (which was The copyright holder for this preprint this version posted February 3, 2020. . https://doi.org /10.1101/19012245 doi: medRxiv preprint To this end, we run multiple simulations with different choices ( 0 1 τ < < ), relative to snail peak ( 0 τ = ). In each case, we estimated mean effect of a given τ -strategy, by averaging over the proper time interval [ ] , 1 τ τ + . Figure 6 shows the results (mean w -values) over 6-year control history for different choices of time shift τ . We find overall effect of τ could be significant between 'best' and 'worst' timing, up to 40% improvement (Y6 over Y1) for trigonometric N and to 50% -for peak N. In both cases, we took high amplitude seasonality ( 0.95 a = ). We observe the optimal timing τ is between one quarter and halfseason for type I, and close to a quarter for type II. Here season start is identified with peak snail population. The results are summarized in Table 2, for each τ we compute ( ) ( )

Optimal timing for molluscicide and integrated strategies
Molluscicide is considered a viable option for control of schistosomiasis 27 . In this section we explore optimal timing for snail control (molluscicide), and integrated strategy (MDA + molluscicide), using logistic snail model (1.9) coupled to MacDonald system (1.7).  is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. not certified by peer review) (which was The copyright holder for this preprint this version posted February 3, 2020. . https://doi.org/10.1101/19012245 doi: medRxiv preprint 0.5 τ = (yellow), close to optimal. In both cases, we took high seasonal variability ( 0.9 a = for type I, and 0.6 a = for type II). Application of repeated molluscicide rapidly brings snail system to a periodic pattern (after 1-2 rounds) for function ( ) N t -total population (unaffected by MDA). But infected snail density ( ) y t approaches zero (due to MDA effect on snail infection).
Next we discuss specifically optimal timing to two types of seasonality. For molluscicide alone the results are shown in Table 3 (see supplement). The optimal timing (maximal MWB reduction) falls near start of K t a reaches its maximum. The overall progress ( The MDA-only strategy looks qualitatively similar to the previous case, prescribed N-function ( Figure 6).
The optimal timing is close to mid-season .5 τ ≈ (Table 4 and supplement) is near mid-season. The mean burden reduction is much more significant now, varies in the range 80-97%, compared to molluscicide.
The combined (MDA+ molluscicide) strategy was implemented at fixed snail removal rate, In the first experiment we run simultaneous MDA and molluscicide with different seasonal timing τ . We also looked at then effect of seasonal amplitude a, strong or weak seasonality. The results are shown in Table 5. Overall, the effect of seasonal timing τ is less significant now, compared to molluscicide or MDA alone, though it increases with higher MDA efficacy. We believe this is due to different optimal values of τ for MDA (typical 0.5 type II). The reason for improved control outcomes at higher a, is the reduced sustainability (endemicity) in highly variable environments, as explained above.
. CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. not certified by peer review) (which was The copyright holder for this preprint this version posted February 3, 2020. . https://doi.org/10.1101/19012245 doi: medRxiv preprint In the final experiment we took independent optimal τ , suggested by earlier simulations, 0.5 τ = for MDA, and 0 τ = for molluscicide. The results are shown in Table 6 (upper half), and compared results to the previous case, i.e. simultaneous MDA+ molluscicide schedule Table 5.
Overall, we achieve further improvement (by factor 2-3) compared to best selection of Table 5 (lower half of Table 6). Depending of seasonal amplitude/type, and MDA coverage-efficacy M ε , the reduction can go as low as 0.5% of the baseline MWB.
We also find that adding molluscicide to MDA can bring a significant improvement, when both controls are properly timed.

Conclusions and discussion
The Basic reproduction number 0 R is a dimensionless parameter widely used in infection modeling. It applies to model reduction (via non-dimensionalization), analysis of equilibria, endemic states and control. In simple cases it can be estimated directly from model inputs (transmission, recovery/ loss rates), and provides full information on endemicity or elimination. One can extend it for large systems with multiple compartments, like SWB, but its scope is more limited. In the SWB models, analysis of endemicity/ elimination requires at least 2 dimensionless parameters (see 9  is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. not certified by peer review) (which was The copyright holder for this preprint this version posted February 3, 2020. . https://doi.org/10.1101/19012245 doi: medRxiv preprint control-induced reproductive number would predict elimination, while the system is locked in a limit cycle. One basic conclusion of our analysis is the need to employ direct dynamic simulations in all cases of unsteady (seasonal) transmission.
Such direct simulation allows one to explore, among other optimal timing for MDA, snail control or integrated strategies in variable environment.
For MDA alone we found such optimal timing is close to half-season, which corresponds to minimal carrying capacity for snails. It's also close to minimal snail density, which could be used as more practical measure for optimal MDA timing. We also found its progress (relative burden reduction), over shortterm 6-year control program, could vary markedly, and 'optimal' timing can achieve up to 50% improvement compared to suboptimal' timing.
For molluscicide control, the optimal timing is shifted closer to start of season, i.e. snail peak, but overall effect of molluscicide on burden reduction is less significant than MDA, even at high killing efficacy.
The combined strategy MDA and molluscicide, can give an additional reduction of worm burden, compared to MDA alone. When both interventions are implemented simultaneously, intra-seasonal variability is less significant, since optimal choices for each one differs. But each intervention implemented at its own optimal timing can achieve much better results.
Overall, snail control can add significantly to the program outcome; in some cases it can lead to nearelimination after short (6-year) duration High seasonal variability gives better control outcomes in all cases, as it makes transmission and endemicity less sustainable.
. CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. not certified by peer review) (which was The copyright holder for this preprint this version posted February 3, 2020.         . CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. not certified by peer review) (which was The copyright holder for this preprint this version posted February 3, 2020. . CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. not certified by peer review) (which was The copyright holder for this preprint this version posted February 3, 2020. . https://doi.org/10.1101/19012245 doi: medRxiv preprint

right). Thick red curve marks the boundary between endemic region (above), and zero-infection (below).
. CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. not certified by peer review) (which was The copyright holder for this preprint this version posted February 3, 2020. . https://doi.org/10.1101/19012245 doi: medRxiv preprint . CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. not certified by peer review) (which was The copyright holder for this preprint this version posted February 3, 2020.  . CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. not certified by peer review) (which was The copyright holder for this preprint this version posted February 3, 2020.  . CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. not certified by peer review) (which was The copyright holder for this preprint this version posted February 3, 2020. . https://doi.org/10.1101/19012245 doi: medRxiv preprint

Supplement SWB system
The SWB approach was developed and refined in several papers (8)(9)(10)(22)(23)(24)(25). It has many advantages compared to MacDonald MWB, simple or 'mated' (NB) version. SWB system carries detailed information on burden distribution and egg release by host population, but it requires no prior assumptions on burden distribution, like NB. Within-host worm biology (mating, mortality, fecundity) is naturally accommodated in SWB. While the number of SWB variables can be large (depending on stratification), it has the same parameters as MacDonald. In fact, the MWB variable ( ) In fact, more detailed information on SWB equilibria can be derived from two transmission confidents A, B rather than their product (like 0 R ). So, communities with identical human infection, but different snail environment (hence different pairs A, B) can produce vastly divergent control outcomes (Supplement).
Parameter 0 R of coupled SWB-snail system Here we derive 0 R formula for coupled SWB-snail system, assuming logistic snail population growth, and SEI snail infection dynamics, and nonlinear snail force of infection (FOI) ( 29 ).
The relevant input parameters are listed in Table 1. . CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. not certified by peer review) (which was The copyright holder for this preprint this version posted February 3, 2020. Formula (0.1) could be further extended to demographically structured SWB communities (e.g. children + SAC + adult et al).
Formula (0.1) can be derived from the Jacobian analysis of the coupled system at zero equilibrium.
Specifically, condition 0 1 R = marks a transition from stable to unstable zero (i.e. stable endemic state).
But there is no direct link between 0 R and endemic state, or possible MDA control outcomes.
. CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. not certified by peer review) (which was The copyright holder for this preprint this version posted February 3, 2020.

Modeling setup
Relaxation patterns to endemic equilibria for reduced and complete MacDonald system  . CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. not certified by peer review) (which was The copyright holder for this preprint this version posted February 3, 2020. . CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. not certified by peer review) (which was The copyright holder for this preprint this version posted February 3, 2020. is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. not certified by peer review) (which was The copyright holder for this preprint this version posted February 3, 2020.

MDA control in stationary environments
Comparison of periodic MDA and its average-model prediction for simple MacDonald The results ( Figure E) show elimination path (yellow curve) for the mean-MDA system, predicted by 1 ef R < , while the exact MDA-response curve (blue) was locked in a periodic (limit) cycle pattern, even above critical frequency ( c T T < ).
. CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. not certified by peer review) (which was The copyright holder for this preprint this version posted February 3, 2020.
Parameter space analysis for mated NB MacDonald MacDonald system with mating and NB worm burden requires 2 dimensionless parameters, instead of single 0 R Figure F: (a) (A, B)   is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. not certified by peer review) (which was The copyright holder for this preprint this version posted February 3, 2020. . https://doi.org/10.1101/19012245 doi: medRxiv preprint Figure F(a) shows that 'endemic persistence -elimination' regions in the (A,B) -parameter space-deviate markedly from isocontours of 0 R . This deviation renders it useless for analysis/prediction of long-term outcomes.

Periodic MDA for MacDonald systems
Here we shall use stationary transmission environment, but vary transmission coefficients.
Specifically, we fixed 0 3 R = , and took 4 A-values in the range 2-3.5. Figure G  prevalence falls below the breakpoint. But higher A (large snail-to-human ratio) have stronger rebound and require more MDA rounds. Figure G (b) show MDA histories for the reduced MacDonald system.
Next we compared MDA predictions of the reduced model with the full host-vector MacDonald system, panels (c-f) of Figure G. The two models have comparable predictions at low transmission 2 A = , but they diverge markedly as A increases. In most cases the reduced model exhibits stronger post-MDA rebound, and therefore takes longer to eliminate, or to reach a breakpoint. The extreme case is

A =
(panel (f)), where complete model goes to elimination after 25 years of MDA, while the reduced model is locked in a limit cycle . CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. not certified by peer review) (which was The copyright holder for this preprint this version posted February 3, 2020. . CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. not certified by peer review) (which was The copyright holder for this preprint this version posted February 3, 2020. is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. not certified by peer review) (which was The copyright holder for this preprint this version posted February 3, 2020. . https://doi.org/10.1101/19012245 doi: medRxiv preprint  Table 1, and amplitude range 0 .9 α ≤ ≤ (blue to red).
. CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. not certified by peer review) (which was The copyright holder for this preprint this version posted February 3, 2020. . https://doi.org/10.1101/19012245 doi: medRxiv preprint . CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. not certified by peer review) (which was The copyright holder for this preprint this version posted February 3, 2020. . https://doi.org/10.1101/19012245 doi: medRxiv preprint