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

^{1}Department of Biology, McGill University, Montréal, QC, Canada^{2}Center for Discrete Mathematics and Theoretical Computer Science, Princeton, NJ, United States^{3}Institute of Parasitology, McGill University, Sainte-Anne-de-Bellevue, QC, Canada^{4}Department of Mathematical and Statistical Sciences, University of Alberta, Edmonton, AB, Canada

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 (${{R}}_{0}$), 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 (${{R}}_{0}$), 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 $\eta \left(\frac{X}{I}\right)={e}^{-\xi \frac{X}{I}\text{}}$, 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 $\alpha \left(S+\eta \left(\frac{X}{I}\right)I+R\right)\left(1-\frac{S+I+R}{K}\right)$. 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 $\frac{X}{I}$. 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 $\delta \left(\frac{X}{I}\right)={E}\frac{X}{I}$ where ${E}$ 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 $\frac{X}{I}$ (i.e., recovery function $\propto \frac{1}{X/I}$), implying that the recovery rate function is proportional to $\frac{Y}{X/I}$ (i.e., the recovery rate function is $\lambda \frac{Y}{X/I}$ 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 ${\int}_{0}^{\tau}g(\psi )d\psi $. Thus the probability that the fish's immunity is not lost τ time units after recovering from infection is ${\int}_{\tau}^{\infty}g(\psi )d\psi .\text{}$ We assume that at a time *t* − τ, $\lambda \frac{Y(t-\tau )I(t-\tau )}{X(t-\tau )}I(t-\tau )$ 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 $g(\tau ):=\frac{{c}^{n}{\tau}^{n-1}{e}^{-c\tau}}{(n-1)!}$, *n* = 1, 2, 3, …, *c* > 0, and $\frac{n}{c}$ 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 $f\left(\frac{X}{I}\right)$ = $\frac{\frac{\theta X}{I}}{\kappa +\frac{X}{I}}$, 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: $\left(1-\frac{X}{P(Y)I}\right)$, 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 $\omega +d+\delta \left(\frac{X}{I}\right)$. 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.

### 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*) ≥ $0,\text{}\overline{S}(t)\ge \text{}0\text{}$for $S(0)\ge \text{}0,\text{}I(0)\ge \text{}0,\text{}R(0)\ge \text{}0,\text{}Y\text{}(0)\ge \text{}0,\text{}X(0)\ge \text{}0\text{}and\text{}\overline{S}(0)\ge \text{}0$. Note that $\frac{X}{I}=\text{}0,\frac{I}{X}=\text{}0$ and $\frac{Y}{X}=\text{}0$ 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 $\dot{\overline{S}}\text{}\left(\overline{S}=0\right)=\frac{cY}{X}{I}^{2}$. Since *I*(*t*) ≥ 0 and $\frac{X}{Y}\ge \text{}0$ for *t* ≥ 0, *I*(0) ≥ 0 and $\frac{X(0)}{Y(0)}\text{}\ge \text{}0$, we have that $\dot{\overline{S}\text{}}(\overline{S}=0)\ge 0,\text{}t\text{}\ge \text{}0$ and $\overline{S}(0)\ge 0$. The non-negativity of *R*(*t*) for *t* ≥ 0 follows from the integral representation (1) and non-negativity of $\frac{Y(t)}{X(t)}$. From the first equation, we have that

From the non-negativity of $\frac{X}{I},I,R,$ we have that if *S*(0) ≥ 0 implies that *S*(*t*) ≥ 0 for *t* ≥ 0.

*Gyrodactylus*-Free Equilibrium Point (GFE) and ${{R}}_{0}$

For the GFE we have that *X* = 0, implying that *I, R, Y* and $\overline{S}$ are null. Plugging this in System (3), we have that $S=\frac{k(\alpha -d)}{\alpha}$. Since *S* > 0, we have that the GFE exist iff α > *d* and is given by $\left(\frac{k(\alpha -d)}{\alpha},\text{}0,\text{}0,0,0,0\right)$. The linearized system corresponding to this equilibrium point is:

The corresponding eigenvalues are:

λ_{1}, λ_{2}, λ_{3} *and λ*_{5} are all less than zero and λ_{4} is less than zero iff μ < *d* + ω. This leads to the definition of the *Gyrodactylus* basic reproduction number, ${{R}}_{0}:=\frac{\mu}{d+\omega}$. Observe that GFE is locally asymptotically stable iff μ < *d* + ω, i.e., GFE is locally asymptotically stable iff ${{R}}_{0}<1$ and unstable if ${{R}}_{0}$ > 1. Therefore, if ${{R}}_{0}$ < 1, the parasite dies out and if ${{R}}_{0}>1$, the parasite will invade the guppy population. ${{R}}_{0}=1$, is a threshold below which the *Gyrodactylus* dies out and above which there is an outbreak. ${{R}}_{0}$ 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.

**Table 1**. Estimated parameter values used to train and test the model and initial values from experiments.

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 ${{R}}_{0}=\text{}2.63$. Since ${{R}}_{0}>1$ 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 ${{R}}_{0}<1$ and for ${{R}}_{0}>1$ 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 ${{R}}_{0}=0.9(<1)$ (A), ${{R}}_{0}=1.1(>1)$ (B) and ${{R}}_{0}=2.63(>1)$ (C) respectively. When ${{R}}_{0}<1$, the system will stabilize to its *Gyrodactylus* free equilibrium $\left(\frac{k(\alpha -d)}{\alpha},\text{}0,\text{}0,0,0,0\right)$ (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 ${{R}}_{0}>1$, 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**. *Gyrodactylus* dynamics with two different basic reproduction numbers (${{R}}_{0}$). **(A)** ${{R}}_{0}=$ 0.90 and **(B)** ${{R}}_{0}=$ 1.1, and **(C)** ${{R}}_{0}=$ 2.59, demonstrating that parasites will die out with an ${{R}}_{0}$ of less than one, and persist if ${{R}}_{0}$ is greater than one.

### 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**. Comparison of model predictions (solid line) of Guppy-*Gyrodactylus* dynamics using the parameters in Table 1 with measured values from averages of laboratory results (diamonds). **(A–D)** we plot the time course dynamics of the number of susceptible guppies, infected guppies, recovered guppies and *Gyrodactylus* per tank, respectively.

### 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**. Comparison of model predictions of Guppy-*Gyrodactylus* dynamics (solid lines) with averages from laboratory data (diamonds) independent from the data used to estimate the model parameters (Table 1). **(A–D)** we have the time course dynamics of the number of susceptible guppies, infected guppies, recovered guppies and *Gyrodactylus* per tank, respectively.

### 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**. Long-term *Gyrodactylus* dynamics when guppy birth is included in the system. The Figure was generated using the parameter values in Table 1 with α = 0.4/fish/day.

### 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 ${{R}}_{0}$, which indicates that an outbreak will occur when the ${{R}}_{0}\text{}$is greater than one and that the system will stabilize to a *Gyrodactylus-*free equilibrium when ${{R}}_{0}$ is less than one. With our parameters, ${{R}}_{0}\text{}$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.

## Data Availability

All datasets generated for this study are included in the manuscript and/or the Supplementary Files.

## Ethics Statement

The animal study was reviewed and approved by McGill University Animal Care Committee (AUP 2014-7547) in compliance with the Canadian Council on Animal Care.

## Author Contributions

The mathematical model was developed at a meeting between CT, JK, GF, and HW. CT collected the data (either experimentally or from published literature) used for parameter estimates and model fitting, while JK validated the model and conducted the sensitivity analysis. CT wrote the initial draft of the manuscript, with editing input from all co-authors.

## Funding

We acknowledge support from FQRNT Equipe grant (MS and GF), NSERC Discovery grants (MS, GF, and HW), NSERC Postdoctoral Fellowship (JK) and NSERC Alexander Graham Bell Canada Graduate Scholarship (CT) and the McGill Department of Biology which provided travel funding for meetings between co-authors.

## Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

The reviewer FD declared a past co-authorship with one of the authors GF to the handling editor.

## Acknowledgments

We thank DIMACS for providing space to conduct the analyses (partially enabled through support from the National Science Foundation under grant #CCF-1445755) as well as Simon Levin's Lab-Princeton University for the same purpose.

## Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fevo.2019.00307/full#supplementary-material

## References

Anderson, R. M., and May, R. M. (1979). Population biology of infectious diseases: part I. *Nature* 280, 361–367. doi: 10.1038/280361a0

Bakke, T. A., Cable, J., and Harris, P. D. (2007). “The biology of gyrodactylid monogeneans: The russian-doll killers,” in *Advances in Parasitology*, eds R. M. J. R. Baker and D. Rollinson. (Cambridge, MA: Academic Press) 459–460. doi: 10.1016/S0065-308X(06)64003-7

Bakke, T. A., Harris, P. D., and Cable, J. (2002). Host specificity dynamics: observations on gyrodactylid monogeneans. *Int. J. Parasitol.* 32, 281–308. doi: 10.1016/S0020-7519(01)00331-9

Barson, N. J., Cable, J., and Van Oosterhout, C. (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. doi: 10.1111/j.1420-9101.2008.01675.x

Cable, J., and van Oosterhout, C. (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. doi: 10.1111/j.1095-8312.2006.00755.x

Chaves, S. S., Gargiullo, P., Zhang, J. X., Civen, R., Guris, D., Mascola, L., et al. (2007). Loss of vaccine-induced immunity to varicella over time. *N. Engl. J. Med.* 356, 1121–1129. doi: 10.1056/NEJMoa064040

Chitnis, N., Hyman, J. M., and Cushing, J. M. (2008). Determining important parameters in the spread of malaria through the sensitivity analysis of a mathematical model. *Bull. Math. Biol.* 70. 1272–96. doi: 10.1007/s11538-008-9299-0

Cornell, S. J., Isham, V. S., and Grenfell, B. T. (2004). Stochastic and spatial dynamics of nematode parasites in farmed ruminants. *Proc. R. Soc. London B Biol. Sci.* 271, 1243–1250. doi: 10.1098/rspb.2004.2744

Croft, D., Arrowsmith, B., Bielby, J., Skinner, K., White, E., Couzin, I., et al. (2003). Mechanisms underlying shoal composition in the Trinidadian guppy, Poecilia reticulata. *Oikos* 100, 429–438. doi: 10.1034/j.1600-0706.2003.12023.x

Dargent, F., Rolshausen, G., Hendry, A. P., Scott, M. E., and Fussmann, G. F. (2016). Parting ways: parasite release in nature leads to sex-specific evolution of defence. *J. Evolut. Biol.* 29, 23–24. doi: 10.1111/jeb.12758

Dargent, F., Scott, M. E., Hendry, A. P., and Fussmann, G. 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. doi: 10.1098/rspb.2013.2371

Denholm, S. J., Hoyle, A. S., Shinn, A. P., Paladini, G., Taylor, N. G., and Norman, R. 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 ONE* 11:e0169168. doi: 10.1371/journal.pone.0169168

Fraser, B. A., and Neff, B. D. (2010). Parasite mediated homogenizing selection at the MHC in guppies. *Genetica* 138, 273–278. doi: 10.1007/s10709-009-9402-y

Fraser, B. A., Ramnarine, I. W., and Neff, B. D. (2009). Selection at the MHC class IIB locus across guppy (*Poecilia reticulata*) populations. *Heredity* 104, 155–167. doi: 10.1038/hdy.2009.99

Gog, J. R., Pellis, L., Wood, J. L., McLean, A. R., Arinaminpathy, N., and Lloyd-Smith, J. O. (2015). Seven challenges in modeling pathogen dynamics within-host and across scales. *Epidemics* 10, 45–48. doi: 10.1016/j.epidem.2014.09.009

Gordon, S. P., Reznick, D. N., Kinnison, M. T., Bryant, M. J., Weese, D. J., Räsänen, K., et al. (2009). Adaptive changes in life history and survival following a new guppy introduction. *Am. Nat.* 174, 34–45. doi: 10.1086/599300

Grenfell, B., and Harwood, J. (1997). (Meta)population dynamics of infectious diseases. *Trends Ecol. Evol.* 12, 395–399. doi: 10.1016/S0169-5347(97)01174-9

Grether, G. F., Kasahara, S., Kolluru, G. R., and Cooper, E. 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. doi: 10.1098/rspb.2003.2526

Hagenaars, T. J., Donnelly, C. A., and Ferguson, N. M. (2004). Spatial heterogeneity and the persistence of infectious diseases. *J. Theoret. Biol.* 229, 349–359. doi: 10.1016/j.jtbi.2004.04.002

Heffernan, J. M., and Keeling, M. J. (2009). Implications of vaccination and waning immunity. *Proc. R. Soc. B Biol. Sci*. 76, 2071–2080. doi: 10.1098/rspb.2009.0057

Heins, D. C., Baker, J. A., Toups, M. A., and Birden, L. S. (2010). Evolutionary significance of fecundity reduction in threespine stickleback infected by the diphyllobothriidean cestode *Schistocephalus solidus*. *Biol. J. Linn. Soc*. 100, 835–846. doi: 10.1111/j.1095-8312.2010.01486.x

Johnsen, B. 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.

Johnsen, B. O., and Jenser, A. J. (1991). The Gyrodactylus story in Norway. *Aquaculture* 98, 289–302. doi: 10.1016/0044-8486(91)90393-L

Kaiser, H., Paulet, T. G., and Endemann, F. (1998). Water quality fluctuations in a closed recirculating system for the intensive culture of the guppy, Poecilia reticulata (Peters). *Aquacul. Res.* 29, 611–615. doi: 10.1046/j.1365-2109.1998.00257.x

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

Kolluru, G. R., Grether, G. F., South, S. H., Dunlop, E., Cardinali, A., Liu, L., 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. doi: 10.1111/j.1095-8312.2006.00675.x

Kuang, Y. (1993). *Delay Differential Equations: With Applications in Population Dynamics*. Arizona: Academic Press.

Lester, R. J. (1972). Attachment of Gyrodactylus to Gasterosteus and host response. *J. Parasitol.* 58, 717–722. doi: 10.2307/3278299

May, R. M., and Anderson, R. M. (1979). Population biology of infectious diseases: part II. *Nature* 280, 455–461. doi: 10.1038/280455a0

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

Mossong, J., Nokes, D. J., Edmunds, W. J., Cox, M. J., Ratnam, S., and Muller, C. P. (1999). Modeling the impact of subclinical measles transmission in vaccinated populations with waning immunity. *Am. J. Epidemiol.* 150, 1238–1249. doi: 10.1093/oxfordjournals.aje.a009951

Okosun, K. O., Ouifki, R., and Marcus, N. (2011). Optimal control analysis of a malaria disease transmission model that includes treatment and vaccination with waning immunity. *Biosystems* 106, 136–145. doi: 10.1016/j.biosystems.2011.07.006

Perez-Jvostov, F., Hendry, A. P., Fussmann, G. F., and Scott, M. E. (2012). Are host–parasite interactions influenced by adaptation to predators? A test with guppies and Gyrodactylus in experimental stream channels. *Oecologia* 170, 77–88. doi: 10.1007/s00442-012-2289-9

Pérez-Jvostov, F., Hendry, A. P., Fussmann, G. F., and Scott, M. E. (2015). Testing for local host–parasite adaptation: an experiment with Gyrodactylus ectoparasites and guppy hosts. *Int. J. Parasitol.* 45, 409–417. doi: 10.1016/j.ijpara.2015.01.010

Reznick, D., and Endler, J. A. (1982). The impact of predation on life history evolution in Trinidadian Guppies (*Poecilia reticulata*). *Evolution* 36, 160–177. doi: 10.1111/j.1558-5646.1982.tb05021.x

Richards, E. L., van Oosterhout, C., and Cable, J. (2010). Sex-specific differences in shoaling affect parasite transmission in guppies. *PLoS ONE* 5:e13285. doi: 10.1371/journal.pone.0013285

Richards, E. L., van Oosterhout, C., and Cable, J. (2012). Interactions between males guppies facilitates the transmission of the monogenean ectoparasite *Gyrodactylus turnbulli*. *Exp. Parasitol.* 132, 483–486. doi: 10.1016/j.exppara.2012.09.019

Richards, G. R., and Chubb, J. 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. doi: 10.1007/s004360050481

Robertson, S., Bradley, J. E., and MacColl, A. 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. doi: 10.1016/j.fsi.2016.11.058

Rosà, R., Pugliese, A., Villani, A., and Rizzoli, A. (2003). Individual-based vs. deterministic models for macroparasites: host cycles and extinction. *Theoret. Popul. Biol.* 63, 295–307. doi: 10.1016/S0040-5809(03)00021-2

Rose, S. M. (1959). Population control in guppies. *Am. Midland Natural.* 62, 474–481. doi: 10.2307/2422539

Sato, A., Figueroa, F., O'Huigin, C., Reznick, D. N., and Klein, J. (1995). Identification of major histocompatibility complex genes in the guppy, P*oecilia reticulata*. *Immunogenetics* 43, 38–49. doi: 10.1007/BF00186602

Scott, M. E. (1982). Reproductive potential of *Gyrodactylus bullatarudis* (Monogenea) on guppies (*Poecilia reticulata*). *Parasitology* 85, 217–236. doi: 10.1017/S0031182000055207

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

Scott, M. 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 D. Rollinson and R. M. Anderson. (New York, NY: Academic Press), 21–38.

Scott, M. E., and Anderson, R. M. (1984). The population dynamics of *Gyrodactylus bullatarudis* (Monogenea) within laboratory populations of the fish host *Poecilia reticulata*. *Parasitology* 89, 159–194. doi: 10.1017/S0031182000001207

Scott, M. E., and Robinson, M. A. (1984). Challenge infections of *Gyrodactylus bullatarudis* (Monogenea) on guppies, *Poecilia reticulata* (Peters), following treatment. *J. Fish Biol.* 24, 581–586. doi: 10.1111/j.1095-8649.1984.tb04828.x

Tadiri, C. P., Dargent, F., and Scott, M. E. (2013). Relative host body condition and food availability influence epidemic dynamics: a *Poecilia reticulata-Gyrodactylus turnbulli* host-parasite model. *Parasitology* 140, 1–9. doi: 10.1017/S0031182012001667

Tadiri, C. P., Scott, M. E., and Fussmann, G. F. (2016). Impact of host sex and group composition on parasite dynamics in experimental populations. *Parasitology* 143:523. doi: 10.1017/S0031182016000172

Tadiri, C. P., Scott, M. E., and Fussmann, G. F. (2018). Microparasite dispersal in metapopulations: a boon or bane to the host population? *Proc. R. Soc. B Biol. Sci.* 285:1885. doi: 10.1098/rspb.2018.1519

van Boven, M., de Melker, H. E., Schellekens, J. F., and Kretzschmar, M. (2000). Waning immunity and sub-clinical infection in an epidemic model: implications for pertussis in The Netherlands. *Math. Biosci.* 164, 161–182. doi: 10.1016/S0025-5564(00)00009-2

Van Oosterhout, C., Harris, P. D., and Cable, J. (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. doi: 10.1046/j.1095-8312.2003.00203.x

Van Oosterhout, C., Mohammed, R. S., Hansen, H., Archard, G. A., McMullan, M., Weese, D. J., et al. (2007). Selection by parasites in spate conditions in wild Trinidadian guppies (*Poecilia reticulata*). *Int. J. Parasitol.* 37, 805–812. doi: 10.1016/j.ijpara.2006.12.016

van Oosterhout, C., Potter, R., Wright, H., and Cable, J. (2008). Gyro-scope: An individual-based computer model to forecast gyrodactylid infections on fish hosts. *Int. J. Parasitol.* 38, 541–548. doi: 10.1016/j.ijpara.2007.09.016

Keywords: epidemic dynamics, mathematical model, guppy, *Gyrodactylus*, host-parasite interactions

Citation: Tadiri CP, Kong JD, Fussmann GF, Scott ME and Wang H (2019) A Data-Validated Host-Parasite Model for Infectious Disease Outbreaks. *Front. Ecol. Evol.* 7:307. doi: 10.3389/fevo.2019.00307

Received: 22 May 2019; Accepted: 30 July 2019;

Published: 21 August 2019.

Edited by:

Pavel Kindlmann, Charles University, CzechiaReviewed by:

Cock Van Oosterhout, University of East Anglia, United KingdomFelipe Dargent, University of Ottawa, Canada

Copyright © 2019 Tadiri, Kong, Fussmann, Scott and Wang. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

*Correspondence: Christina P. Tadiri, christina.tadiri@mail.mcgill.ca; Jude D. Kong, jk1567@dimacs.rutgers.edu; Hao Wang, hao8@ualberta.ca

^{†}These authors have contributed equally to this work as first authors