A Data-Validated Host-Parasite Model for Infectious Disease Outbreaks

Abstract

The use of model experimental systems and mathematical models is important to further understanding of infectious disease dynamics and strategize disease mitigation. Gyrodactylids are helminth ectoparasites of teleost fish which have many dynamical characteristics of microparasites but offer the advantage that they can be quantified and tracked over time, allowing further insight into within-host and epidemic dynamics. In this paper, we design a model to describe host-parasite dynamics of the well-studied guppy-Gyrodactylus turnbulli system, using experimental data to estimate parameters and validate it. We estimate the basic reproduction number (), for this system. Sensitivity analysis reveals that parasite growth rate, and the rate at which the guppy mounts an immune response have the greatest impact on outbreak peak and timing both for initial outbreaks and on longer time scales. These findings highlight guppy population average resistance and parasite growth rate as key factors in disease control, and future work should focus on incorporating heterogeneity in host resistance into disease models and extrapolating to other host-parasite systems.

Introduction

The guppy-Gyrodactylus system is a well-known model host-parasite system, used in numerous experimental and field studies (Scott, 1985a,b; Richards and Chubb, 1998; Cable and van Oosterhout, 2007). Guppies, Poecilia reticulata, are a common ovoviviparous tropical teleost fish whose abundance and ability to survive a broad range of environmental variables and availability in pet stores worldwide have made them ideal subjects for research in various disciplines. Gyrodactylus spp. (Monogenea) are ectoparasites which feed on the epithelial cells and mucus of many marine and freshwater teleost fish species (Bakke et al., 2007). They attach to the epidermis of their host via specialized hooks and are directly transmitted primarily by jumping to a new host during contact (Scott and Anderson, 1984; Kearn, 1994). They also reproduce directly on the host, with the developing embryo containing within itself a second developing embryo, which allows for rapid population growth of the parasite directly on an infected host (Kearn, 1994; Bakke et al., 2007). Upon infection, hosts mount an immune response, including mucus secretion (Lester, 1972), as well as a non-specific complement which kills gyrodactylids (Sato et al., 1995; Woo, 2006; Bakke et al., 2007; Robertson et al., 2017). Gyrodactylid infection can result in high rates of mortality (Van Oosterhout et al., 2003), and induce a temporary refractory period in surviving hosts (Scott and Robinson, 1984; Scott, 1985a,b).

In general, parasites are typically divided into two categories: microparasites (such as viruses and bacteria) which are microscopic and tend to proliferate and transmit rapidly, often leading to high morbidity and mortality and inducing acquired resistance in surviving hosts and therefore causing periodic epidemics, while macroparasites (such as worms and insects) often have more complex life cycles and persist in populations, often with overdispersed distributions among hosts and rarely causing severe morbidity and mortality or acquired resistance (Anderson and May, 1979; May and Anderson, 1979). Due to their rapid growth rate and infection-induced refractory period, gyrodactylids cause periodic epidemic outbreaks, making their population dynamics typical of microparasites like viruses and bacteria (Anderson and May, 1979) despite being helminths which traditionally fall into the macroparasite category. However, they have a key distinction from other typical microparasites in that parasite population size or burden is a central factor in determining host-parasite dynamics, as it directly influences transmission, mortality and several other parameters. Intrinsic population dynamics of Gyrodactylus sp. on isolated fish have been identified under standardized environmental conditions (Scott, 1982, 1985a,b). Both short- and long-term dynamics of Gyrodactylus sp. within laboratory populations of guppies have also been observed (Scott, 1985a,b; Richards and Chubb, 1998; Tadiri et al., 2013, 2016). However, the need for a comprehensive model that can describe and make predictions for this system and others like it still exists.

Traditional microparasite Susceptible, Infected, Recovered (SIR) models can effectively describe epidemic dynamics of infectious diseases for which the parasite population size is unknown or less relevant than host category of infection (Anderson and May, 1979; Grenfell and Harwood, 1997; Hagenaars et al., 2004). Yet, SIR models are less applicable to parasites such as Gyrodactylus spp., where parasite burden plays a crucial role in host-parasite population dynamics. Although macroparasite models directly consider parasite number, they also often overlook dynamics in parasite numbers within individual hosts (Anderson and May, 1979; Rosà et al., 2003; Cornell et al., 2004). As gyrodactylids and many other parasites do not fit neatly into the micro-/ macro-parasite dichotomy, there is a clear need for a unifying framework which considers both host and parasite populations (Gog et al., 2015). Previous efforts to mathematically describe this system using various types of models have captured basic initial epidemic dynamics, but failed to effectively describe longer-term fluctuations due to gradual loss of immunity over time (Scott and Anderson, 1984; van Oosterhout et al., 2008). Similarly, infection dynamics on individual fish have been simulated, but the broader scale transmission and population dynamics were not incorporated (van Oosterhout et al., 2008). The objective for this paper is to establish a mathematical model that effectively describes experimental data on guppy-Gyrodactylus dynamics, particularly with regards to host immunity waning and longer-term dynamics and to estimate the sensitivity to various parameters that we have not been able to effectively test in the laboratory. Ideally, this model can be applied to other directly transmitted, directly reproducing parasites for which parasite burden impacts host-parasite relations, with a waning immunity post-infection.

Methods

The guppy-Gyrodactylus system shares some key characteristics with common directly transmitted infectious disease dynamics (e.g., infection-induced host mortality, refractory period, infection by host-to-host contact). Therefore, we design an SIR-type model with distributed delay (which captures the varying immunity period of the guppy) to describe the dynamics of guppies and Gyrodactylus. Since the guppy immune response plays a crucial role in eliminating Gyrodactylus, we explicitly integrate the dynamics of the immune response into the model. Thereafter, the distributed delay model is converted to an equivalent system of ordinary differential equations using the linear chain approach. Next, we ensure that non-negative initial values do not give rise to a negative solution. To determine the Gyrodactylus basic reproduction number (), the stability analysis of the Gyrodactylus-free equilibrium point was performed. This threshold is particularly of use because it allows us to determine the maximum potential number of Gyrodactylus that will be produced due to the introduction of one Gyrodactylus in a Gyrodactylus-free population of guppies, which can help inform control measures. Model parameters unavailable in the literature were estimated by data fitting using previously published experimental data. Next, the model was validated by comparison to measurements from independent but analogous laboratory experiments. Finally, using the estimated parameters, together with the parameters from the literature, the sensitivity of the outbreak peak magnitude and the time to outbreak peak to the parameters of the model was determined. This sensitivity could be useful in determining the most influential parameters for designing control measures.

Derivation of the Model

In this section we derive a guppy-Gyrodactylus interaction model with distributed delay. Figure 1 provides a conceptual flowchart for the system. The model consists of four coupled equations tracing the rates of change of guppy population (G), guppy immune response (Y) and Gyrodactylus population (X). The guppy total population is divided into three sub-groups: susceptible (S), infected (I), and recovered (R) guppies. The change in number of susceptible guppies could be due to (1) birth by any guppy (we assume all guppies are born susceptible), (2) loss of immunity by a recovered guppy, (3) death of a susceptible guppy or (4) becoming infected due to contact with an infected guppy. We assume the guppy population to be homogenous and the natural birth rate of the guppy is assumed to be constant, α. The birth rate of infected individuals is diminished by a function that is linearly proportional to the number of parasite it harbors, which we assumed to be , where ξ is the steepness of parasite-induced fecundity reduction. Although it is unclear whether guppy fecundity is reduced by Gyrodactylus infection, this is the case for many infectious diseases, including those of fish (Heins et al., 2010) so we allow for it in our model, while defaulting the parameters to 0 in our simulations because in our experimental populations no birth was observed. Instead of the exponential growth assumed in typical SIR-type models, we consider a logistic growth for the guppy population because uninfected guppies exhibit density-dependent population growth up to a carrying capacity K (Rose, 1959). Thus, the growth rate of the guppy is . The rate at which susceptible fish become infected, is described by the function:

where b is a constant. The natural death rate of guppies is assumed to be a constant d. The parasite mean intensity is represented as . The population of infected guppies can increase when an infected guppy contacts a susceptible guppy, resulting in transmission, and decreases when any one of them dies or recovers. In addition to the natural death rate of guppies, infected guppies may also be killed by Gyrodactylus at a rate described by the function where is a constant. The recovery rate function is assumed to be directly proportional to the average immune response Y (i.e., recovery function ∝ Y) and inversely proportional to the parasite intensity (i.e., recovery function ), implying that the recovery rate function is proportional to (i.e., the recovery rate function is where λ is the proportionality constant). Immunity is assumed to affect both the rate of parasite population growth and the rate at which infected fish become recovered and recovered fish regain susceptibility. Every recovered guppy is assumed to acquire an immunity that wanes with time following initial infection. To model immunity waning, we assume that the immunity period of every fish varies from 0 to ∞ in order to capture the wide variability in the period of acquired resistance observed among many guppy populations (Scott and Robinson, 1984; Scott, 1985a,b; Richards and Chubb, 1998; Cable and van Oosterhout, 2007). We let g(ψ) denote the probability density function that a fish takes exactly ψ time units to lose its immunity after recovering from infection, which implies that the probability that a fish's immunity is lost τ time units after recovering from infection, is . Thus the probability that the fish's immunity is not lost τ time units after recovering from infection is We assume that at a time t − τ, fish left the infected population compartment and joined the recovered population. The probability that these fish are still alive τ time units after leaving the infected compartment is e^{−d τ}. Hence the total number of recovered fish at time t, is:

We consider g(τ) to be the density function for a gamma distribution , n = 1, 2, 3, …, c > 0, and to be the average duration of the immune memory. As the number of parasites increases, the guppy immune system gradually builds a defense against the parasite at a rate = , proportional to the density of non-specific immune complement responsible for killing Gyrodactylus, where θ is the maximum rate of increase of immunity. This defense gradually reduces the guppy parasite carrying capacity. The per capita growth rate of the parasite follows a logistic growth: , where p(Y) is the average parasite carrying capacity of an infected guppy. p(Y) is assumed to be an exponentially decreasing function of the average immune response of the guppy p(Y) = re^−γY, where r and γ are constants. The natural parasite death rate is assumed to be a constant, ω. We assume that when a guppy dies, all the parasites on it die, since dead guppies are more likely to be predated or washed downstream (Van Oosterhout et al., 2007). In our experiments, dead fish were removed from tanks < 1 day after death in order to minimize transmission from dead fish, but it's possible some could have occurred. Thus, the total per capita parasite death rate is . Thus, we have the following systems of delay differential equations:

where ϕ_i, i = 1, 2, 3, 4 and 5 are bounded continuous functions, μ is the maximum per capita parasite growth rate. We assume that.

Figure 1

Reduction to an Ordinary Differential Equation Model

We assume that g(τ) = ce^−cτ (i.e. n = 1) and apply the chain trick method (Kuang, 1993) to convert System (2) to a system of ordinary differential equations: Let

Substituting this in System (2) reduces it to the following system of ODE:

Positivity and Basic Reproduction Number

In this section, we show that non-negative initial data give rise to non-negative solutions, establish conditions for the existence and stability of the Gyrodactylus-free equilibrium point of the system and determine the basic reproduction number.

Positivity

Positivity and boundedness of a model guarantee that the model is biologically well-behaved. For positivity of the System (3), we have the following theorem:

Theorem 1. All solutions of System (3) are positive for all t in (0, ∞)

Proof. We need to show that S(t) ≥ 0, I(t) ≥ 0, R(t) ≥ 0, Y (t) ≥ 0, X(t) ≥ for . Note that and for X < 1 or I < 1. We start by proving that if I(0) ≥ 0, ⇒ I(t) ≥ 0 for all t > 0. From the second equation of Systems (3), we see that İ(I = 0) = 0. Thus I(t) ≥ 0 for t ≥ 0. Also, from the third equation, we have that Ẋ(X = 0) = 0 for X(0) ≥ 0. Hence X(t) ≥ 0 for t ≥ 0. From the last equation, we have that . Since I(t) ≥ 0 and for t ≥ 0, I(0) ≥ 0 and , we have that and . The non-negativity of R(t) for t ≥ 0 follows from the integral representation (1) and non-negativity of . From the first equation, we have that

From the non-negativity of we have that if S(0) ≥ 0 implies that S(t) ≥ 0 for t ≥ 0.

Gyrodactylus-Free Equilibrium Point (GFE) and

For the GFE we have that X = 0, implying that I, R, Y and are null. Plugging this in System (3), we have that . Since S > 0, we have that the GFE exist iff α > d and is given by . The linearized system corresponding to this equilibrium point is:

The corresponding eigenvalues are:

λ₁, λ₂, λ₃and λ₅ are all less than zero and λ₄ is less than zero iff μ < d + ω. This leads to the definition of the Gyrodactylus basic reproduction number, . Observe that GFE is locally asymptotically stable iff μ < d + ω, i.e., GFE is locally asymptotically stable iff and unstable if > 1. Therefore, if < 1, the parasite dies out and if , the parasite will invade the guppy population. , is a threshold below which the Gyrodactylus dies out and above which there is an outbreak. has an intuitive biological interpretation: it is the average number of Gyrodactylus resulting from the introduction of a single Gyrodactylus into an otherwise Gyrodactylus-free population over the course of its life span.

Parameter Estimation and Model Validation Using Independent Measurements

We used data obtained from separate laboratory previously published experiments to, respectively, estimate the model parameters not available in the literature and to test the fit of the model. To estimate the model parameters we used experimental data averaged from four groups of eight male fish where one fish per group was infected with two parasites and the infection was allowed to spread naturally throughout the tank (Tadiri et al., 2018). These fish were bred in the lab from pet-store “feeder” guppies. To test the model, we used experimental data averaged from four groups of eight fish (four males and four females) where parasites were introduced to each group via a donor juvenile fish infected with three parasites that was removed once at least three parasites had naturally transferred to the experimental fish (Tadiri et al., 2016). These fish were third-generation lab reared fish bred from 33 originally family lines originally obtained from wild population in Trinidad but mixed haphazardly in experimental tanks. In both cases, 2-3 parasites were introduced in order to keep the introduction as close to one as possible while minimizing the probability of accidental parasite death or that an old or male parasite would be introduced preventing reproduction. No difference in host-parasite dynamics was found among all-male, all-female and mixed sex groups of eight fish (Tadiri et al., 2016).

In both experiments, each fish was individually marked, and the number of parasites on each fish was counted every other day to obtain SIR numbers and total parasite population size. In our all experimental groups, no birth was observed within the 42 days. Hence, we ignore vital dynamics for guppies in the model.

To estimate the parameter values of System (3), we use the non-linear regression function nlinfit(.) in MATLAB. The function nlinfit(.) uses the Levenberg-Marquardt algorithm (Moré, 1978) to fit the solution of the biodegradation module to the data. Some of parameters used in solving System (3) namely ω, d, and K, were taken from the literature (Rose, 1959; Scott and Anderson, 1984): the units, values and source of these parameters are provided in Table 1.

The validity of our model in predicting Gyrodactylus outbreak was evaluated by using the estimated parameters in the model to generate S, I, R Gyrodactylus and immune response data then comparing the predicted data to measured data using the goodnessOfFit(.) function in MATLAB.

Sensitivity Analysis

The objective of this subsection is to discuss the sensitivity of the magnitude of the initial Gyrodactylus outbreak peak, and time to the initial outbreak peak to the parameters of the system. For this analysis, we use the normalized forward sensitivity index (Chitnis et al., 2008):

where F^* is the quantity being considered.

Since we do not have the explicit formula for the initial outbreak peak, or time to peak, we use central difference approximation to estimate them:

Letting h = 1% of the parameter value (P), Equation (4) becomes:

Longer-Term Dynamics and Generic Sensitivity Analysis

In this section, we simulate the longer-term parasite dynamics in a system that allows for guppy vital dynamics. We use an average birth rate estimated from literature of 0.4/fish/day (Rose, 1959) with no infection-induced reduction in fecundity. We equally assess the sensitivity of the generic outbreak peaks and period to the parameters of the system from the long-term simulation.

Results

System Basic Reproduction Number and Gyrodactylus-Free Equilibrium Point

Using the parameters in Table 1 we have that . Since for the estimated parameters values, the GFE is not asymptotically stable, meaning that an introduction of one Gyrodactylus into a naïve guppy population would result in an outbreak. Next, we illustrate the dynamics of the system for and for using values very close to one. Rearranging the fourth equation of System (2), we have:

Figure 2 shows the long-term behavior of the Gyrodactylus population for (A), (B) and (C) respectively. When , the system will stabilize to its Gyrodactylus free equilibrium (Panel A). The number of Gyrodactylus, the number of infected guppies, the number of recovered guppies and the guppy immune compliment density tend to zero as t increases. When , there will be a Gyrodactylus outbreak (Panels B and C). These outcomes are robust for large sets of initial values and parameter values.

Figure 2

Fitting the Model to Data

Table 1 contains the value of the parameters obtained from fitting System (2) to the experimental data described above. Figure 3 shows the simulated susceptible, infected, recovered guppies and Gyrodactylus dynamics along with measured data. We obtained a goodness-of-fit statistic (NMSE) value of 0.99. This statistic indicates that the model is able to predict the training data accurately.

Figure 3

Model Evaluation Using NMSE

Using the parameter values in Table 1, with the procedure described above, we assess the validity of our model in predicting guppy-Gyrodactylus dynamics data. Figure 4 shows a comparison between our simulated and measured data. The goodness-of-fit statistics suggests that System (2) with the given parameter values is a good fit for guppy-Gyrodactylus dynamics data (NMSE = 0.70).

Figure 4

Sensitivity of the Magnitude of the Initial Outbreak Peak

The sensitivity indices of the magnitude of the peak of the initial outbreak measure how the magnitude of the peak of the initial outbreak depends on different parameters. Table 2 contains the sensitivity indices of the amplitude of the first outbreak peak obtained using Equation (5). The two parameters with the greatest independent influence on the system were parasite increase rate (μ), and maximum rate of increase of immunity (θ).

Sensitivity of the Time to Initial Outbreak Peak

Sensitivity indices of the time to initial outbreak peak measure how the first epidemic outbreak time depends on different parameters as seen in the Table 3. Similar to the sensitivity of the magnitude of the peak of the first outbreak, the parasite increase rate (μ) and the maximum rate of the increase of immunity (θ) were the most influential parameters.

Longer-Term Dynamics

By allowing for natural guppy birth in our systems, we are able to simulate steady state oscillating dynamics such as those expected of epidemic infectious diseases with generic peaks and periods (Figure 5).

Figure 5

Sensitivity of the Generic Outbreak Peak Magnitude

Sensitivity indices of the generic outbreak peak measure how the magnitude of subsequent outbreak peaks under steady state oscillating dynamics depend on different parameters (Table 4). The guppy immune related parameters, θ (maximum rate of increase of immunity and v (rate of decay of immunity) have the strongest relationship to the outbreak peak. The negative value of the sensitivity of the outbreak to θ, indicates that a low value of θ will lead to a more severe parasite outbreak. The positive value of the sensitive of the outbreak peak to v, on the other hand, tells us that if the guppy immunity wanes faster, there will be a severe outbreak.

Sensitivity of the Generic Outbreak Period

Sensitivity indices of the generic outbreak period measure how the time to subsequent outbreak peaks under steady state oscillating dynamics depend on different parameters (Table 5). The parasite increase rate (μ) and the maximum rate of increase of immunity (θ) have the strongest relationship to the outbreak period.

Discussion

In this paper, we define a mathematical model that effectively describes guppy-Gyrodactylus dynamics in small populations. In estimating parameters based on both literature and our own experimental data we determined our model to accurately describe the dynamics of this system. Additionally, we validated our model using a neutral data set from a separate experiment and found that it fit reasonably well. We also find the model to be mathematically and biologically sound through our analysis of , which indicates that an outbreak will occur when the is greater than one and that the system will stabilize to a Gyrodactylus-free equilibrium when is less than one. With our parameters, was greater than one, indicating that an outbreak will occur in our system with the introduction of one parasite. This model builds of previous efforts to model this system (Scott and Anderson, 1984), but incorporates a more realistic representation of immunity. Firstly, and most importantly, we describe fish immune response to infection in its own equation, rather than assuming a linear constant to represent immunity. We specifically also describe the waning of immunity post-infection using a distributed delay function, which allows for repeated, dampening cycles of outbreaks without constant immigration of naive hosts. We also allow for host population growth rather than fixed immigration and consider a parasite-induced reduction in fecundity (Perez-Jvostov et al., 2012), which had not been investigated at the time of previous models. Longer-term simulations using population growth estimates from literature with our other parameter estimates from experiments demonstrate that our model is capable of describing oscillating parasite dynamics typical of those observed in the wild (Van Oosterhout et al., 2007). These developments are important to more accurately explaining guppy-Gyrodactylus dynamics and could have a broader applicability to other systems as well.

Gyrodactylus are a large genus of over 400 ectoparasites infecting at least 20 orders of teleost fish (Bakke et al., 2002), and our model can most directly be applied to other species of this genus. Gyrodactylids have had significant economic impact, causing epizootics in many resources fish such as carp, trout and African catfish (Woo, 2006) and, most notably Atlantic salmon fisheries, particularly in Norway in the 1960's and ‘70’s which saw large declines due to G. salaris (Johnsen, 1978; Bakke et al., 2007), and efforts to recover these populations and prevent disease spread to other watersheds are still ongoing (Denholm et al., 2016). Since the basic life cycle of gyrodactylids and their relationships with their hosts are similar (Bakke et al., 2002, 2007) this model would only need reparameterization to be applied to a range of other aquaculture species. Beyond other gyrodactylids, many infectious diseases also confer immunity that decays over time. Our methods of applying a distributed delay to describe waning guppy immunity to Gyrodactylus are novel to this system, and can also be used for other infectious diseases with declining immunity, most notably being comparable to the waning of vaccine-induced immunity which has observed in many human diseases (Heffernan and Keeling, 2009) such as measles (Mossong et al., 1999), pertussis (van Boven et al., 2000), malaria (Okosun et al., 2011), and varicella (Chaves et al., 2007) and modeled using different methods. Given the broader applicability of our methods, our results have important implications for disease management, as we identify the most impactful parameters on disease outbreaks, and thus crucial intervention points.

Our sensitivity analysis found that the most influential parameters on both initial outbreak amplitude and time to initial outbreak in our system were parasite increase rate (μ) and maximum increase rate of guppy immunity (θ). These results indicate parasite growth rate and host resistance play the strongest role in the severity and speed of an outbreak, which makes logical sense. The higher parasite growth rate, or lower the immune response, the greater the parasite abundance will be. Given the small population sizes of our experiments, it is possible that host density may have affected the relative importance of these variables compared to the transmission rate or population size. The density of fish in our experiments was higher than wild populations (Croft et al., 2003), but lower than in commercial guppy populations (Kaiser et al., 1998), therefore an average approximation of the different conditions in which Gyrodactylus outbreaks may occur. Given the relatively short timescale of our experiments, we did not observe any impacts of longer-term parameters such as guppy birthrate or natural guppy mortality. However, sensitivity analysis of our generic outbreak magnitudes and periods consistently demonstrate that parasite virulence and parameters relating to guppy immunity have the strongest impact on our system and therefore this result was not an artifact of our experimental design. Both guppy resistance (Van Oosterhout et al., 2003; Dargent et al., 2016), and parasite virulence (Cable and van Oosterhout, 2007) are known to evolve rapidly and vary widely among populations due to different selective pressures and our findings indicate that understanding this heterogeneity is of significance to predicting and controlling disease outbreaks.

One limitation of this model is that it was based on laboratory, rather than field data and large differences in both host mortality and parasite burdens have been observed between the lab and field settings. Gyrodactylids persist in the wild and are observed at typically low burdens are observed, however mark-recapture experiments have suggested that infection can cause severe mortality (Van Oosterhout et al., 2007) typical of epidemics and in aquaculture (Johnsen, 1978; Johnsen and Jenser, 1991), and laboratory (Scott and Anderson, 1984; Van Oosterhout et al., 2003) settings, outbreaks are known to cause severe disease and mortality. Also, in the wild, guppies inhabit streams which are punctuated by “pools” separated by waterfalls, thus creating a network of populations among which unidirectional migration of hosts (and potentially parasites) downstream is possible (Van Oosterhout et al., 2007; Barson et al., 2009), however in our current model we focus only on populations in isolation. Furthermore, as previously mentioned, the population sizes were much smaller than those in a natural setting, and as such, the timescales used to estimate some of the longer-term parameters of our model, such as host birth and immunity waning to full susceptibility may not fully reflect dynamics in the wild (Van Oosterhout et al., 2007). Nevertheless, our estimates obtained from both literature and short-term data show a good fit for our data and could potentially be directly applied to aquaculture settings with only reparameterization specific to the species of interest. Moreover, longer-term simulations with our model show its ability to predict longer-term fluctuating dynamics such as those observed in the wild, indicating the predictive value of this model to more natural settings.

Another limitation is our assumption of homogeneity of hosts, which is not accurate in this system. Guppies are known to exhibit a broad range in both life history traits (Reznick and Endler, 1982; Gordon et al., 2009) and innate resistance to parasites (Fraser et al., 2009; Fraser and Neff, 2010; Dargent et al., 2013), both within and among populations. Additionally, individual guppies may vary in their susceptibility to parasites due to individual characteristics such as size (Cable and van Oosterhout, 2007; Tadiri et al., 2013), carotenoid coloration (Grether et al., 2004; Kolluru et al., 2006) and sex (Richards et al., 2010, 2012; Dargent et al., 2016; Tadiri et al., 2016). Our parameters don't capture the wide variability that occurs in nature, or how this heterogeneity may influence host-parasite dynamics in the population but were instead based on average values obtained from literature and our own laboratory observations. Moreover, significant variability even in average population-level resistance has been observed among wild populations and domestic fish to various strains of gyrodactylids (Van Oosterhout et al., 2003; Cable and van Oosterhout, 2007; Dargent et al., 2013, Pérez-Jvostov et al., 2015) and it's possible that our estimated parameters may not fit some extreme cases of particularly low- or high-resistance population. However, despite not accounting for such complexities, our model fit data from two experiments, one which used fish from various wild populations from Trinidad and one which used domestic fish, therefore we find these average values to be a decent approximation.

In conclusion, we were able to develop and validate a mathematical model that more effectively describes the guppy-Gyrodactylus system, thus contributing to a further understanding of disease dynamics. Through sensitivity analysis, we were able to identify key factors affecting outbreaks to strategize control measures for parasites which increase in numbers due to reproduction directly on the host (in the absence of transmission) and are directly transmitted via host contact, particularly those relating to parasite growth rate and host resistance. Our findings have implications for a broader range of systems, with our model being most directly applicable to other gyrodactylids such as G. salaris, which is known to cause severe mortality and morbidity in Atlantic salmon fisheries (Johnsen, 1978; Bakke et al., 2007), but these methods could be also applicable to many other infections for which immunity decays over time, such as that observed for some vaccines.

References

• 1

  AndersonR. M.MayR. M. (1979). Population biology of infectious diseases: part I. Nature280, 361–367. 10.1038/280361a0

  • CrossRef
  • Google Scholar

• 2

  BakkeT. A.CableJ.HarrisP. D. (2007). The biology of gyrodactylid monogeneans: The russian-doll killers, in Advances in Parasitology, eds BakerR. M. J. R.RollinsonD.. (Cambridge, MA: Academic Press) 459–460. 10.1016/S0065-308X(06)64003-7

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 3

  BakkeT. A.HarrisP. D.CableJ. (2002). Host specificity dynamics: observations on gyrodactylid monogeneans. Int. J. Parasitol.32, 281–308. 10.1016/S0020-7519(01)00331-9

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 4

  BarsonN. J.CableJ.Van OosterhoutC. (2009). Population genetic analysis of microsatellite variation of guppies (Poecilia reticulata) in Trinidad and Tobago: evidence for a dynamic source–sink metapopulation structure, founder events and population bottlenecks. J. Evolut. Biol.22, 485–497. 10.1111/j.1420-9101.2008.01675.x

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 5

  CableJ.van OosterhoutC. (2007). The role of innate and acquired resistance in two natural populations of guppies (Poecilia reticulata) infected with the ectoparasite Gyrodactylus turnbulli. Biol. J. Linn. Soc.90, 647–655. 10.1111/j.1095-8312.2006.00755.x

  • CrossRef
  • Google Scholar

• 6

  ChavesS. S.GargiulloP.ZhangJ. X.CivenR.GurisD.MascolaL.et al. (2007). Loss of vaccine-induced immunity to varicella over time. N. Engl. J. Med.356, 1121–1129. 10.1056/NEJMoa064040

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 7

  ChitnisN.HymanJ. M.CushingJ. M. (2008). Determining important parameters in the spread of malaria through the sensitivity analysis of a mathematical model. Bull. Math. Biol.70. 1272–96. 10.1007/s11538-008-9299-0

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 8

  CornellS. J.IshamV. S.GrenfellB. T. (2004). Stochastic and spatial dynamics of nematode parasites in farmed ruminants. Proc. R. Soc. London B Biol. Sci.271, 1243–1250. 10.1098/rspb.2004.2744

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 9

  CroftD.ArrowsmithB.BielbyJ.SkinnerK.WhiteE.CouzinI.et al. (2003). Mechanisms underlying shoal composition in the Trinidadian guppy, Poecilia reticulata. Oikos100, 429–438. 10.1034/j.1600-0706.2003.12023.x

  • CrossRef
  • Google Scholar

• 10

  DargentF.RolshausenG.HendryA. P.ScottM. E.FussmannG. F. (2016). Parting ways: parasite release in nature leads to sex-specific evolution of defence. J. Evolut. Biol.29, 23–24. 10.1111/jeb.12758

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 11

  DargentF.ScottM. E.HendryA. P.FussmannG. F. (2013). Experimental elimination of parasites in nature leads to the evolution of increased resistance in hosts. Proc. R. Soc. B Biol. Sci.280:2371. 10.1098/rspb.2013.2371

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 12

  DenholmS. J.HoyleA. S.ShinnA. P.PaladiniG.TaylorN. G.NormanR. A. (2016). Predicting the potential for natural recovery of atlantic salmon (Salmo salar L.) Populations following the Introduction of Gyrodactylus salaris Malmberg, 1957 (Monogenea). PLoS ONE11:e0169168. 10.1371/journal.pone.0169168

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 13

  FraserB. A.NeffB. D. (2010). Parasite mediated homogenizing selection at the MHC in guppies. Genetica138, 273–278. 10.1007/s10709-009-9402-y

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 14

  FraserB. A.RamnarineI. W.NeffB. D. (2009). Selection at the MHC class IIB locus across guppy (Poecilia reticulata) populations. Heredity104, 155–167. 10.1038/hdy.2009.99

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 15

  GogJ. R.PellisL.WoodJ. L.McLeanA. R.ArinaminpathyN.Lloyd-SmithJ. O. (2015). Seven challenges in modeling pathogen dynamics within-host and across scales. Epidemics10, 45–48. 10.1016/j.epidem.2014.09.009

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 16

  GordonS. P.ReznickD. N.KinnisonM. T.BryantM. J.WeeseD. J.RäsänenK.et al. (2009). Adaptive changes in life history and survival following a new guppy introduction. Am. Nat.174, 34–45. 10.1086/599300

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 17

  GrenfellB.HarwoodJ. (1997). (Meta)population dynamics of infectious diseases. Trends Ecol. Evol.12, 395–399. 10.1016/S0169-5347(97)01174-9

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 18

  GretherG. F.KasaharaS.KolluruG. R.CooperE. L. (2004). Sex–specific effects of carotenoid intake on the immunological response to allografts in guppies (Poecilia reticulata). Proc. R. Soc. London. Series B Biol. Sci.271, 45–49. 10.1098/rspb.2003.2526

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 19

  HagenaarsT. J.DonnellyC. A.FergusonN. M. (2004). Spatial heterogeneity and the persistence of infectious diseases. J. Theoret. Biol.229, 349–359. 10.1016/j.jtbi.2004.04.002

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 20

  HeffernanJ. M.KeelingM. J. (2009). Implications of vaccination and waning immunity. Proc. R. Soc. B Biol. Sci.76, 2071–2080. 10.1098/rspb.2009.0057

  • CrossRef
  • Google Scholar

• 21

  HeinsD. C.BakerJ. A.ToupsM. A.BirdenL. S. (2010). Evolutionary significance of fecundity reduction in threespine stickleback infected by the diphyllobothriidean cestode Schistocephalus solidus. Biol. J. Linn. Soc.100, 835–846. 10.1111/j.1095-8312.2010.01486.x

  • CrossRef
  • Google Scholar

• 22

  JohnsenB. O. (1978). The effect of an attack by the parasite Gyrodactylus salaris on the population of salmon parr in the river Lakselva, Misvaer in northern Norway. J. Arct. Biol.11, 7–9.

  • Google Scholar

• 23

  JohnsenB. O.JenserA. J. (1991). The Gyrodactylus story in Norway. Aquaculture98, 289–302. 10.1016/0044-8486(91)90393-L

  • CrossRef
  • Google Scholar

• 24

  KaiserH.PauletT. G.EndemannF. (1998). Water quality fluctuations in a closed recirculating system for the intensive culture of the guppy, Poecilia reticulata (Peters). Aquacul. Res.29, 611–615. 10.1046/j.1365-2109.1998.00257.x

  • CrossRef
  • Google Scholar

• 25

  KearnG. C. (1994). Evolutionary expansion of the Monogenea. Int. J. Parasitol.24, 1227–1271. 10.1016/0020-7519(94)90193-7

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 26

  KolluruG. R.GretherG. F.SouthS. H.DunlopE.CardinaliA.LiuL.et al. (2006). The effects of carotenoid and food availability on resistance to a naturally occurring parasite (Gyrodactylus turnbulli) in guppies (Poecilia reticulata). Biol. J. Linn. Soc.89, 301–309. 10.1111/j.1095-8312.2006.00675.x

  • CrossRef
  • Google Scholar

• 27

  KuangY. (1993). Delay Differential Equations: With Applications in Population Dynamics.Arizona: Academic Press.

  • Google Scholar

• 28

  LesterR. J. (1972). Attachment of Gyrodactylus to Gasterosteus and host response. J. Parasitol.58, 717–722. 10.2307/3278299

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 29

  MayR. M.AndersonR. M. (1979). Population biology of infectious diseases: part II. Nature280, 455–461. 10.1038/280455a0

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 30

  MoréJ. J. (1978). The Levenberg-Marquardt algorithm: implementation and theory. Numer. Anal.630, 105–116. 10.1007/BFb0067700

  • CrossRef
  • Google Scholar

• 31

  MossongJ.NokesD. J.EdmundsW. J.CoxM. J.RatnamS.MullerC. P. (1999). Modeling the impact of subclinical measles transmission in vaccinated populations with waning immunity. Am. J. Epidemiol.150, 1238–1249. 10.1093/oxfordjournals.aje.a009951

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 32

  OkosunK. O.OuifkiR.MarcusN. (2011). Optimal control analysis of a malaria disease transmission model that includes treatment and vaccination with waning immunity. Biosystems106, 136–145. 10.1016/j.biosystems.2011.07.006

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 33

  Perez-JvostovF.HendryA. P.FussmannG. F.ScottM. E. (2012). Are host–parasite interactions influenced by adaptation to predators? A test with guppies and Gyrodactylus in experimental stream channels. Oecologia170, 77–88. 10.1007/s00442-012-2289-9

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 34

  Pérez-JvostovF.HendryA. P.FussmannG. F.ScottM. E. (2015). Testing for local host–parasite adaptation: an experiment with Gyrodactylus ectoparasites and guppy hosts. Int. J. Parasitol.45, 409–417. 10.1016/j.ijpara.2015.01.010

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 35

  ReznickD.EndlerJ. A. (1982). The impact of predation on life history evolution in Trinidadian Guppies (Poecilia reticulata). Evolution36, 160–177. 10.1111/j.1558-5646.1982.tb05021.x

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 36

  RichardsE. L.van OosterhoutC.CableJ. (2010). Sex-specific differences in shoaling affect parasite transmission in guppies. PLoS ONE5:e13285. 10.1371/journal.pone.0013285

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 37

  RichardsE. L.van OosterhoutC.CableJ. (2012). Interactions between males guppies facilitates the transmission of the monogenean ectoparasite Gyrodactylus turnbulli. Exp. Parasitol.132, 483–486. 10.1016/j.exppara.2012.09.019

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 38

  RichardsG. R.ChubbJ. C. (1998). Longer-term population dynamics of Gyrodactylus bullatarudis and G. turnbulli (Monogenea) on adult guppies Poecilia reticulata in 50-l experimental arenas. Parasitol. Res.84, 753–756. 10.1007/s004360050481

  • CrossRef
  • Google Scholar

• 39

  RobertsonS.BradleyJ. E.MacCollA. D. C. (2017). No evidence of local adaptation of immune responses to Gyrodactylus in three-spined stickleback (Gasterosteus aculeatus). Fish Shellf. Immunol.60, 275–281. 10.1016/j.fsi.2016.11.058

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 40

  RosàR.PuglieseA.VillaniA.RizzoliA. (2003). Individual-based vs. deterministic models for macroparasites: host cycles and extinction. Theoret. Popul. Biol.63, 295–307. 10.1016/S0040-5809(03)00021-2

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 41

  RoseS. M. (1959). Population control in guppies. Am. Midland Natural.62, 474–481. 10.2307/2422539

  • CrossRef
  • Google Scholar

• 42

  SatoA.FigueroaF.O'HuiginC.ReznickD. N.KleinJ. (1995). Identification of major histocompatibility complex genes in the guppy, Poecilia reticulata. Immunogenetics43, 38–49. 10.1007/BF00186602

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 43

  ScottM. E. (1982). Reproductive potential of Gyrodactylus bullatarudis (Monogenea) on guppies (Poecilia reticulata). Parasitology85, 217–236. 10.1017/S0031182000055207

  • CrossRef
  • Google Scholar

• 44

  ScottM. E. (1985a). Dynamics of challenge infections of Gyrodactylus bullatarudis Turnbull (Monogenea) on guppies, Poecilia reticulata (Peters). J. Fish Dis.8, 495–503. 10.1111/j.1365-2761.1985.tb00964.x

  • CrossRef
  • Google Scholar

• 45

  ScottM. E. (1985b). Experimental epidemiology of Gyrodactylus bullatarudis (Monogenea) on guppies (Poecilia reticuata: short- and long-term studies, in Ecology and Genetics of Host-Parasite Interactions, eds RollinsonD.AndersonR. M.. (New York, NY: Academic Press), 21–38.

  • Google Scholar

• 46

  ScottM. E.AndersonR. M. (1984). The population dynamics of Gyrodactylus bullatarudis (Monogenea) within laboratory populations of the fish host Poecilia reticulata. Parasitology89, 159–194. 10.1017/S0031182000001207

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 47

  ScottM. E.RobinsonM. A. (1984). Challenge infections of Gyrodactylus bullatarudis (Monogenea) on guppies, Poecilia reticulata (Peters), following treatment. J. Fish Biol.24, 581–586. 10.1111/j.1095-8649.1984.tb04828.x

  • CrossRef
  • Google Scholar

• 48

  TadiriC. P.DargentF.ScottM. E. (2013). Relative host body condition and food availability influence epidemic dynamics: a Poecilia reticulata-Gyrodactylus turnbulli host-parasite model. Parasitology140, 1–9. 10.1017/S0031182012001667

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 49

  TadiriC. P.ScottM. E.FussmannG. F. (2016). Impact of host sex and group composition on parasite dynamics in experimental populations. Parasitology143:523. 10.1017/S0031182016000172

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 50

  TadiriC. P.ScottM. E.FussmannG. F. (2018). Microparasite dispersal in metapopulations: a boon or bane to the host population?Proc. R. Soc. B Biol. Sci.285:1885. 10.1098/rspb.2018.1519

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 51

  van BovenM.de MelkerH. E.SchellekensJ. F.KretzschmarM. (2000). Waning immunity and sub-clinical infection in an epidemic model: implications for pertussis in The Netherlands. Math. Biosci.164, 161–182. 10.1016/S0025-5564(00)00009-2

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 52

  Van OosterhoutC.HarrisP. D.CableJ. (2003). Marked variation in parasite resistance between two wild populations of the Trinidadian guppy, Poecilia reticulata (Pisces: Poeciliidae). Biol. J. Linnean Soc.79, 645–651. 10.1046/j.1095-8312.2003.00203.x

  • CrossRef
  • Google Scholar

• 53

  Van OosterhoutC.MohammedR. S.HansenH.ArchardG. A.McMullanM.WeeseD. J.et al. (2007). Selection by parasites in spate conditions in wild Trinidadian guppies (Poecilia reticulata). Int. J. Parasitol.37, 805–812. 10.1016/j.ijpara.2006.12.016

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 54

  van OosterhoutC.PotterR.WrightH.CableJ. (2008). Gyro-scope: An individual-based computer model to forecast gyrodactylid infections on fish hosts. Int. J. Parasitol.38, 541–548. 10.1016/j.ijpara.2007.09.016

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 55

  WooP. T. (2006). Fish Diseases and Disorders.Oxfordshire, UK: CABI Publishing.

  • Google Scholar