Experimental investigation of alternative transmission functions: Quantitative evidence for the importance of nonlinear transmission dynamics in host–parasite systems

Abstract Understanding pathogen transmission is crucial for predicting and managing disease. Nonetheless, experimental comparisons of alternative functional forms of transmission remain rare, and those experiments that are conducted are often not designed to test the full range of possible forms. To differentiate among 10 candidate transmission functions, we used a novel experimental design in which we independently varied four factors—duration of exposure, numbers of parasites, numbers of hosts and parasite density—in laboratory infection experiments. We used interactions between amphibian hosts and trematode parasites as a model system and all candidate models incorporated parasite depletion. An additional manipulation involving anaesthesia addressed the effects of host behaviour on transmission form. Across all experiments, nonlinear transmission forms involving either a power law or a negative binomial function were the best‐fitting models and consistently outperformed the linear density‐dependent and density‐independent functions. By testing previously published data for two other host–macroparasite systems, we also found support for the same nonlinear transmission forms. Although manipulations of parasite density are common in transmission studies, the comprehensive set of variables tested in our experiments revealed that variation in density alone was least likely to differentiate among competing transmission functions. Across host–pathogen systems, nonlinear functions may often more accurately represent transmission dynamics and thus provide more realistic predictions for infection.

Density-dependent transmission (typically the default option for most models of disease; McCallum et al., 2001) assumes rates of host-to-host contact increase linearly with population density, whereas frequency-dependent (density-independent) transmission (typically used to model sexually transmitted-and vector-borne infections; McCallum et al., 2001) assumes a constant per-capita rate of contact with other hosts regardless of population density (Begon et al., 2002). These models make strikingly different predictions for infection dynamics (Hoch, Fourichon, Viet, & Seegers, 2008;Smith et al., 2009;Wonham, Lewis, Rencławowicz, & van den Driessche, 2006), in particular regarding whether there is a threshold host population size for pathogen persistence (Bolker & Greenfell, 1995;Swinton, Harwood, Grenfell, & Gilligan, 1998); density-dependent transmission predicts such a threshold, whereas frequency-dependent transmission does not, with the consequence that parasites showing frequency-dependent transmission can in theory drive the host population extinct (De Castro & Bolker, 2005). Given these differences, making accurate predictions about how diseases affect host population dynamics depends on identifying the appropriate transmission function, with important implications for disease management (Jennelle et al., 2014;McCallum et al., 2001;Smith et al., 2009).
Despite the ubiquity of density-or frequency-dependent transmission functions in disease ecology theory, it is clear that, though mathematically convenient, these functions do not adequately capture the functional form of transmission seen in nature, and in some cases, their underlying assumptions have been proven invalid (Cross et al., 2013;Morters et al., 2013;Smith et al., 2009;Wonham et al., 2006). Given this, a range of alternative transmission functions has been proposed (McCallum et al., 2001; Table 1). Many of these functions assume that at low density contacts are likely to be directly proportional to host and parasite density (i.e. approximating density-dependent transmission), but saturate at very high host or parasite densities (approximating frequency-dependent transmission; Antonovics, Iwasa, & Hassell, 1995;D'Amico, Elkinton, Dwyer, Burand, & Buonaccorsi, 1996;Fenton et al., 2002;McCallum et al., 2001). Such a pattern has previously been found from analysis of infection incidence data of cowpox virus in wild rodents (Smith et al., 2009), finding that transmission was best described by a hybrid function which moved between density-dependent transmission when host densities were low, and frequency-dependent (density-independent) transmission when host densities were high. A similar hybrid function has also been found to provide a better description of elk contact rates and the spread of brucellosis than either frequency or density dependence (Cross et al., 2013).
Rather than representing categorical alternatives, density and frequency dependence may capture two specific points on a wider spectrum of possible transmission shapes (McCallum et al., 2017).
Reflecting this, more flexible, phenomenological functions have been proposed (e.g. power functions) that allow transmission to take a range of nonlinear forms beyond the density-/frequency-dependent extremes (Fenton et al., 2002;Hochberg, 1989Hochberg, , 1991Liu, Levin, & Iwasa, 1986; Table 1). Mechanisms that can give rise to such nonlinear transmission functions include heterogeneities in the distribution of infectious particles, or density-dependent mortality of the pathogen (Briggs & Godfray, 1995). However, beyond a general sense that transmission is likely to be more nonlinear than assumed by standard formulations, a comprehensive understanding of how this is driven by variation in host and parasite density or abundance is lacking (McCallum et al., 2017).
Empirically it is challenging to quantify the overall magnitude of transmission, let alone determine its functional form. Several studies have sought to do this using infection data from natural systems by embedding alternative transmission functions within larger population dynamic models and then fitting those models to observed incidence data (Begon et al., 1998(Begon et al., , 1999Cross et al., 2013;Rachowicz & Briggs, 2007). However, those approaches may be subject to uncertainty arising from potential inaccuracies in the specification of the wider population dynamic model or quantifying host densities, including the relevant area over which to evaluate transmission and the occurrence or timing of infection. An alternative approach is to conduct transmission experiments in the laboratory, where many of these factors can be controlled or quantified (D'Amico et al., 1996;Knell, Begon, & Thompson, 1996;Ryder, Miller, White, Knell, & Boots, 2007). Typically these experiments involve exposing different densities of uninfected (susceptible) and infected hosts (or parasite infective stages, depending on the mode of transmission), and quantifying the number of hosts that become infected over a period of time. However, the specific experimental design is crucial to differentiating among alternative transmission functions; many functions are indistinguishable in certain contexts, for example if increasing host abundance also increases host density (i.e. if arena size remains constant). Hence, distinguishing among competing transmission functions often requires independent manipulation of both host and parasite numbers and arena size (May & Anderson, 1979;McCallum et al., 2001), which surprisingly few studies have explored.
Here we experimentally tested a suite of mathematical functions as competing hypotheses to represent the transmission dynamics of a model system involving interactions between amphibian hosts and trematode parasites. Macroparasite transmission (i.e. arthropods and helminths) involves processes analogous to microparasite transmission (i.e. bacteria and viruses) with added mechanistic control often well suited for experiments. Because macroparasites often infect hosts through free-living stages, each of which represent independent and quantifiable infection events, this system allows us to decouple the host and parasite components of transmission to help understand transmission in systems with these characteristics. Specifically, we used a maximum likelihood approach to determine which mathematical transmission function was best supported by experiments that systematically varied duration of exposure, host density, parasite density and the total number of parasites. Furthermore, we examined the role of host behaviour (via experimental manipulation using anaesthesia). Finally, we tested the generality of our conclusions by analysing two other datasets for trematode transmission (Karvonen, Paukku, Valtonen, & Hudson, 2003;Paller, Kimura, & Uga, 2007), ultimately finding that nonlinear forms of transmission were more appropriate functions than either the classic density-dependent functions, or density-independent and ratio-dependent functions analogous to frequency-dependent transmission for macroparasites.

| Study system
Ribeiroia ondatrae is a trematode in the family Echinostomatidae with a complex life cycle, sequentially infecting snails (Helisoma trivolvis), larval amphibians and finally amphibian-eating birds (Johnson, Sutherland, Kinsella, & Lunde, 2004;Tkach, Kudlai, & Kostadinova, 2016). Transmission of R. ondatrae from snails to amphibians occurs through direct infection by free-living aquatic parasite stages-cercariae-which then form encysted metacercariae in the amphibian that induce developmental malformations (Johnson, Lunde, Ritchie, & Launer, 1999). During the transmission process, cercariae exhibit a searching behaviour until they locate a suitable encystment site (Beaver, 1939;Johnson et al., 2004;Sutherland, 2005). Preferred encystment locations in amphibians include just below the epidermis and above the muscular tissues at the base of hindlimbs, inguinal and tail resorption area, cloaca and T A B L E 1 The suite of transmission functions used to model transmission of the free-living infective stages (cercariae, C) to encysted stages (metacercariae, M) of Ribeiroia ondatrae in Pacific chorus frog (Pseudacris regilla) tadpoles. Transmission is defined successful acquisition of parasite infective stages to the host over time in units of numbers of metacercariae (M). The form of each function used to model microparasite transmission from the literature is provided with the form used for macroparasites in this study for comparison. In the microparasite functions, S is the number of susceptible individuals (analogous to hosts, H) and I is the number of infectious individuals (analogous to cercariae, C). In all functions, β is the transmission parameter, assumed here to be constant in time. Additionally, v is the volume of the enclosure, p and q are the susceptible (host) and infectious (cercariae) responses that represent how densities of each independently affect transmission efficiency. Finally, k is the time-dependent index of aggregation parameter for the negative binomial model. References refer to microparasite functions lower mandible (Beaver, 1939;Johnson & Hartson, 2009;Johnson et al., 2004;Sutherland, 2005). We used the Pacific chorus frog, Pseudacris regilla, in the family Hylidae as our focal amphibian host because it is among the species with the highest reported frequencies of malformations, approaching 90% of emerging metamorphs in some populations (Johnson, Preston, Hoverman, & Richgels, 2013). Mortality and pathology in the amphibian host is related to the intensity of infection-the total number of parasites in the host-highlighting the need for understanding transmission dynamics (Johnson et al., 2012).

| Model development and evaluation
We compared 10 candidate models representing unique hypotheses for pathogen transmission with experimental data using an information theoretic approach following the methods of Greer, Briggs, and Collins (2008) and Rachowicz and Briggs (2007). We defined transmission in terms of successful acquisition of parasite infective stages to the host over time. We developed our set of models (Table 1) Table 1. For example, ϕ = βΗ represents the constant risk function for hosts, where H is the number of hosts and β is the transmission parameter encompassing contact between C and H and the probability of infection. We ignored other sources of cercariae mortality (Karvonen et al., 2003;Paller et al., 2007), but these were likely negligible over the short time frame of our experiments (15-240 min). We represent populations of hosts and parasites as numbers with density included explicitly, when applicable, by incorporating the volume (v) of the experimental tanks into the functional forms (Begon et al., 2002).
For most of the transmission functions listed in Table 1, we found the analytical solutions to each differential equation. For a subset (negative binomial 1 and 2), there were no analytical solutions, so we solved them numerically using the deSolve package (Soetaert, Petzoldt, & Setzer, 2010) in r (R Core Team, 2015). There are two different formulations for the negative binomial model used in the literature, one with hosts and one without (Briggs & Godfray,1995;Greer et al., 2008;Rachowicz & Briggs, 2007

| Laboratory transmission experiments
We used a series of targeted, laboratory experiments that manipulated different variables independently with the aim of testing the alternative transmission functions ( and (4)  Finally, our experimental design minimized heterogeneity in hosts and parasites by using tadpoles that were size-and stage-matched and cercariae collected within a narrow age range (Appendices S1 and S2), thereby limiting the mechanisms that could drive nonlinear transmission dynamics.

| Host behaviour experiment
To examine the role of host behaviour in influencing the form of transmission, we performed a fifth experiment in which we reduced host activity by anaesthesia with neutral-buffered MS-222 (0.125% tricaine methanesulfonate for 3 min; Appendix S3). By anaesthetizing tadpoles with a vertebrate neurotoxin, we eliminated anti-parasite behaviours such as evasive manoeuvres, while avoiding a direct influence on trematode cercariae as supported by our observations and by previous studies (Daly & Johnson, 2011;Koprivnikar, Forbes, & Baker, 2006). Furthermore, previous research found no effect of MS-222 on tadpole resistance to trematode infection (Sears, Snyder, & Rohr, 2013). We followed the same procedures used for varying parasite number while maintaining parasite density across container volumes (see above), with four to seven replicates per container (12 tadpoles recovered early from anaesthesia and were excluded from analysis). After each experimental exposure, we maintained tadpoles for 48 hr before quantifying parasite infection (Appendix S5). We also evaluated the effect of anaesthesia on the probability of cercarial encystment using an overdispersed binomial generalized linear model (Crawley, 2007;R Core Team, 2015).

| Evaluation of other host-macroparasite systems
We tested the generality of our findings with R. ondatrae and amphibian hosts by extending our analysis to previously published data from two additional macroparasite systems: rainbow trout (Oncorhynchus mykiss) and the trematode Diplostomum spathaceum (Karvonen et al., 2003) and minnow (Zacco temmincki) and the trematode Centrocestus armatus (Paller et al., 2007). These studies experimentally exposed individual hosts to different (1) cercariae densities and (2) cercariae numbers independent of density (Appendix S5), but without any manipulation of host density or exposure duration. Although the original studies evaluated frequency-and density-dependent functions, our reanalysis fit all models described in Table 1.

| Overview of laboratory transmission experiments
The manipulations produced a wide range of R. ondatrae metacercariae per host that scaled with variation in parasite, host or time treatments. Duration of exposure and host density had especially strong effects on transmission success (Figure 1a,b). Increases in exposure duration had a strong positive effect leading to the maximum infection level overall with c. 80% of cercariae successfully encysting after 240 min. Increasing host density led to an overall greater transmission success for cercariae, but individual hosts had lower average infection (metacercariae per host) when multiple tadpoles were present.
With respect to differentiating among the candidate transmission functions, variation in the total number of parasites (independent of density [cercariae/L]), variation in host density and variation in the duration of exposure offered the clearest differences, while variation in parasite density alone discriminated poorly ( Figure 1, Table 2). Across most experiments, we found strong support for two types of nonlinear functions, namely the power law and negative binomial functions, and poor support for classical density-dependent, density-independent and ratio-dependent models (Table 2).

| Varying parasite number (constant parasite density, variable volume)
Increasing the number of parasites while keeping parasite density constant led to a nonlinear saturation of infection ( Figure 1a). For instance, average infection per tadpole increased from 2.8 ± 1.1 (SD) metacercariae when exposed to four cercariae to 15.7 ± 9.2 (SD) when exposed to 144 cercariae. This relationship was best represented by the negative binomial function, and both negative binomial models were equally supported because they reduced to the same model with fixed H (Figure 2a, Table 2a). Overall, there was a sharp decline in the proportion of successful cercariae from 0.7 to 0.1 over the range of volumes tested (0.12-4.8 L).

| Varying host density (constant volume, variable host number)
Increases in host density led to a slightly nonlinear increase in total cercariae transmission success (Figure 1b), which increased from 37% (23.1 ± 7.3 recovered parasites of an initial 63) with a single host to 51% (32.1 ± 7.7) with four tadpoles. This was the largest in- Models without a host component received little support (Table 2).

| Varying duration of exposure
Increasing the duration of exposure caused a rapid increase in the proportion of successful cercariae from 0.3 to over 0.7 followed by saturation of infection at ~0.8 over exposure periods longer than 60 min (Figure 1c).
This corresponded to an average infection of 8.3 ± 2.7 (SD) metacercariae per tadpole after 15 min rising to 21.8 ± 3.6 (SD) at 60 min. Based on this nonlinear saturation, the power law C model showed the greatest support, with the power law CH model next best with ΔAIC c = 2.3 (Table 2c). All remaining functions underestimated the rapid transmission dynamics initially and overestimated the total infection at saturation (ΔAIC c >183).

| Varying parasite density (constant volume, variable parasite number)
Finally, we observed a linear increase in average infection with increasing parasite density from 1.4 ± 1.2 (SD) (28% of total parasites) at an exposure of five cercariae to 13.8 ± 11.0 (SD) with exposure to 60 cercariae (23% of total parasites) (Figure 1d). This experiment was least able to differentiate among the competing transmission functions, providing almost equal support across all the tested functions (Table 2d). The proportion of successful cercariae remained relatively constant over the range of parasite densities from ~0.17 to 0.28.

| Evaluation of other host-macroparasite systems
Our analysis of data from previous studies of macroparasite transmission showed congruence with our empirical results, suggesting that nonlinear transmission functions may be general across a variety of parasite and host taxa and different scales of experimental procedures (Table S1). Importantly, the results for both systems contrast the original conclusions that transmission was frequency-dependent (density-independent) (Figure 3; Karvonen et al., 2003;Paller et al., 2007). However, Karvonen et al. (2003) log-transformed the data prior to analysis, which may have led to the support for a linear function for otherwise nonlinear dynamics. Similar to our empirical data, experiments manipulating parasite density were less able to distinguish between functional forms compared to when parasite numbers were varied independent of density.
For rainbow trout (O. mykiss) and D. spathaceum cercariae (Karvonen et al., 2003), there was a rapid increase in transmission followed by saturation at the highest exposure levels when the number of parasites increased independently of parasite density. This relationship supported the power C model, unlike in our experiment where the results conformed more closely with the negative binomial model ( Figure 3a, Table S1). For the experiment where parasite density was varied, the negative binomial and power law C models were almost equally supported (ΔAIC c 0.3-3); however, the trends appeared nearly linear similar to our experimental results from varying parasite density ( Figure 3b). For minnows (Z. temmincki) and C. armatus cercariae (Paller et al., 2007), the negative binomial model was most strongly supported when parasite numbers were varied independent of parasite density, and when parasite density varied (Figure 3c,d). This function captured the apparent curvilinear relationship when parasite numbers increased and aligned closely with results of the current study involving R. ondatrae and amphibian hosts.

| DISCUSSION
By using an array of experimental manipulations that independently varied host density, parasite density, total parasite exposure and exposure duration, we tested the efficacy of 10 different functional forms of host-parasite transmission that broadly encompass previous research from both epidemiological and ecological literatures (Table 1). We found that nonlinear transmission forms involving either a power law or negative binomial function were the best-fitting models, consistently outperforming both density-dependent and density-independent functions. Thus, transmission of the trematode R. ondatrae to amphibian hosts within controlled exposure trials scaled in a positive but nonlinear fashion with the number of infectious parasites, the number of hosts and the duration of exposure.
With the exception of the parasite density experiment where all models were nearly indistinguishable, all formulations of densitydependent and density-independent transmission were very poor Nonlinear transmission forms related to host and parasite densities received strong support in our analysis. Past applications of nonlinear transmission include modelling of insect pathogens (Briggs & Godfray, 1995), chytrid infection in natural populations of Daphnia (Johnson, Ives, Lathrop, & Carpenter, 2009) and Schistosoma sp. cercariae (Gao, Liu, Luo, & Xie, 2011). Mechanisms proposed to help explain nonlinear transmission dynamics include spatial aggregation in the distribution of hosts or parasites, heterogeneity in host susceptibility (immunity and physiology), and variation in the behaviour of hosts or parasitoids (Briggs & Godfray, 1995;Hochberg, 1991 (Beaver, 1939;Haas, 1994). Moreover, potential interference by large numbers of simultaneously infecting cercariae may contribute to the F I G U R E 3 The number of Diplostomum spathaceum metacercariae infecting rainbow trout (Oncorhynchus mykiss) across a range of (a) parasite numbers (C), and (b) parasite density (/L) (C/v), (data from Karvonen et al., 2003). The number of Centrocestus armatus metacercariae attached to minnow (Zacco temmincki) across a range of (c) parasite numbers (C), and (d) parasite density (/L) (C/v), (data from Paller et al., 2007). Points represent the infection levels of individual fish. Lines represent expected values for the best-fitting models (Table S1). For some functions and experimental conditions, models reduced to the same functional form (a, d) resulting in exactly the same Akaike information criterion (AIC These two perspectives are not the same, and together they combine to determine the shape of the overall transmission function, further emphasizing the need to independently vary both host and parasite density and number to capture the full interaction between them. Despite the importance of accurately modelling transmission for forecasting disease, there are relatively few empirical studies of transmission for macroparasites, limiting opportunities to develop a generalized framework. Parasites with free-living infective stages are often assumed to have the same transmission mode as directly transmitted pathogens, or they are not explicitly modelled because of the short time-scale of the dynamics (May & Anderson, 1979;Rachowicz & Briggs, 2007). Omitting particular experimental conditions or only including classical forms often leads to support for frequencydependent transmission even when reanalysis indicated that nonlinear functions were a better fit (Karvonen et al., 2003;Paller et al., 2007).
For instance, in a recent reanalysis of data involving interactions between the human blood fluke Schistosoma mansoni and its snail intermediate hosts, Civitello and Rohr (2014) found that mechanistic models of transmission produced nonlinear dynamics, and power law models were superior to traditional density-dependent transmission, which did not account for parasite depletion.
Nonlinear functions lead to complex dynamics in terms of stability and the thresholds for pathogen invasion in host populations (Hochberg, 1991;Liu et al., 1986), and suggest there may also be a threshold density of infected hosts or infective stages required for disease persistence (Knell et al., 1996). Unsuitable transmission functions can therefore lead to inaccurate conclusions about the risk of pathogen-induced extinction, the effectiveness of control strategies and the periodicity of host-pathogen cycles (De Castro & Bolker, 2005;Greer et al., 2008;Hochberg, 1989;Morters et al., 2013;Smith et al., 2009). Because of their flexibility, nonlinear forms of transmission may be more useful across a variety of host-pathogen systems in forecasting disease occurrence and number of hosts infected as opposed to strictly frequency-or density-dependent transmission or when the true functional form for transmission is unknown (Fenton et al., 2002). Extending the current models to larger spatial scales and longer temporal scales will facilitate exploration of questions related to host populations and communities, including the role of changes in biodiversity, as well as the efficacy of proposed management strategies (Dobson, 2004;Johnson et al., 2013;Rudolf & Antonovics, 2005).
By considering other ecological interactions represented by nonlinear functions, transmission dynamics can be integrated with concepts from consumer resource dynamics and natural enemy ecology (Lafferty et al., 2015). In our models, the transmission coefficient is analogous to the searching efficiency of parasitoids (Knell et al., 1996). Transmission functions are also analogous to the functional response of predators (McCallum et al., 2001). This suggests that cercariae searching behaviour might be similar to parasitoids or predators (Combes, Bartoli, & Theron, 2002;Fenton & Rands, 2004). One intriguing possibility is the potential interference between cercariae for preferential encystment locations.
At the short time-scales at which transmission occurs, this mechanism may limit infection at the highest exposure levels. Further evidence for this competitive interaction between cercariae could be the increased variance in infection we observed in both the high parasite density and parasite number experiments, which have more potential for parasite interactions with variable outcomes. Our models addressed average infection, but additional experiments with larger sample sizes could help identify mechanisms for the variance and potential for parasite aggregation among hosts .