## Abstract

All animals and plants are, to some extent, susceptible to disease caused by varying combinations of parasites, viruses and bacteria. In this paper, we present a mathematical model of interactions between a host, two parasitoids and a pathogen which shows that the presence of an infection can preserve and promote diversity in such multi-species systems. Initially, we use a system of ordinary differential equations to investigate interactions between two species of parasitoids, a host and a host infection. We show that the presence of all four species is necessary for the system as a whole to persist, and that in particular, the presence of the pathogen is necessary for the coexistence of the two parasitoid species. The inclusion of infection induces a wide range of dynamics, including chaos, and these dynamics are robust for a wide range of parameter values. We then extend the model to include spatial effects by introducing random motility (diffusion) of all three species and examine the subsequent spatio-temporal dynamics, including travelling waves and other more complicated heterogeneous behaviour. The computational simulation results of the model suggest that infection in the hosts can blunt the effects of competition between parasitoids, allowing the weaker competitor to survive. Regardless of the nature of the stability of the coexistent steady state of the system, there is an initial period of transient dynamics, the length of which can be extended by an appropriate choice of initial conditions. The existence of these transient dynamics suggests that systems subject to regular restoration to a starting state, such as agro-ecosystems, may be kept in a continual state of dynamic transience, and this has implications for the use of natural enemies to control insect pests, the preservation of biodiversity in farmland habitats and the more general dynamics of disease processes.

## 1. Introduction

The impact of disease on natural populations is well recognized (Anderson & May 1991; Delfino & Simmons 2000), and simple ecological models have been frequently employed in attempts to assess the likely impact of pathogens on host populations (Anderson & May 1980, 1981). Since then, there have been particularly significant advances in developing both system-specific models (Briggs & Godfray 1995, 1996) and in understanding the natural ecology of host–pathogen interactions (Grenfell & Dobson 1995; Hudson *et al*. 2002; Bonsall 2004). It is also clear that both pathogens and parasites can affect the dynamics of populations, even at low prevalence (Anderson 1995; Tomkins & Begon 1999). When combined with other factors, such as competition or predation, interspecific infection can play a crucial role in determining the structure of an ecological community under circumstances where the infection does not affect all species equally (Bowers & Turner 1997; Tomkins *et al*. 2003). However, the indirect effects of infection on competing species have received considerably less attention.

Parasitoids are insects whose larvae develop by feeding on a single arthropod host. They may be solitary, laying only one egg per host, or gregarious laying anything from a few to several thousand eggs per host. Parasitoids are also classified as idiobionts when they immobilize or kill their hosts at the time of oviposition, or as koinobionts when they permit the host to continue and feed, with the host dying only when the juvenile parasitoids emerge from their eggs. They may be specialists, attacking a single species of host or generalists, able to use a variety of host species. One reason why parasitoids have been so extensively studied is that they provide a very accurate and effective experimental metaphor for many exploiter–resource associations in ecology, such as predator–prey and herbivore–plant interactions (Godfray 1994).

In many associations, both hosts and parasitoids are subject to disease, but while there has been extensive ecological research into the effects of disease on host–parasitoid systems, and some work on host–parasitoid–pathogen systems where only one species of parasitoid is present (Hochberg *et al*. 1990), there has been less work on the modelling of such systems (Dobson & Grenfell 1995; Gulland 1995; Sait *et al*. 2000). Hochberg *et al*. (1990) considered the effects of inteference between parasitoid and pathogen in an infected host and examined the dynamics of the host–parasitoid–pathogen interactions. Disease may be transmitted between hosts in a number of ways—horizontally through ingestion of droppings (frass), contact with detritus from infected hosts, through the ‘dirty needle’ effect (Dye *et al*. 1995) where a parasitoid attacks an infected host and then carries the infection to an uninfected host on its ovipositor or vertically through infection of eggs from an infected mother, as in the case of insect baculoviruses (Sait *et al*. 1994; Briggs *et al*. 1995; Burden *et al*. 2003). Many diseases make a host more susceptible to mortality if it is simultaneously subject to other sources of physiological stress, such as food deprivation, while other infections depress the host's immune system or alter its ability to defend itself (Wertheim *et al*. 2005; Fytrou *et al*. 2006).

We consider here a situation in which hosts are infected both vertically through ovarian transmission and horizontally through contact with infection-bearing host products, such as frass. Such a situation might arise among phytophagous moths in the tropics where larvae ingest frass as they consume leaf material but where frequent rain rapidly cleans vegetation or where the infection is otherwise short lived outside the host. The infection is sublethal, but has an indirect detrimental effect on the parasitoids. The mathematical model which will be presented in §2 simulates the effects of such infections in natural communities, therefore, by introducing disease dynamics into a simple model system of two parasitoids attacking a single host species. We undertake an analysis of the steady states of the model and examine their stability as well as investigating transient dynamics using computational simulations. In order to examine the situation where the various species can move freely and interact with each other within a given spatial domain, the model is then extended to include spatial effects in both hosts and parasitoids through the inclusion of random motility (cf. diffusion) in all species. We then present a numerical analysis of the spatio-temporal dynamics of this system. The results of the analysis of the models have implications not only for understanding the role of infection in determining the structure of ecological communities, but also for the effective application of biological control in agro-ecosystems, the impact such methods might have on the ability of farmland habitats to sustain biodiversity and the general dynamics of disease processes.

## 2. The mathematical model

Initially, we consider a system of two species of gregarious parasitoids, denoted *P*_{1} and *P*_{2}, attacking a single species of host denoted *H*. The model constitutes a system of three ordinary differential equations—one for the host population and one for each species of parasitoid—describing the interactions between the species. We assume that the host population has logistic growth, with an intrinsic rate of growth *r* and a carrying capacity *K*. The parameter *r* can be defined in several ways as is appropriate to the system under study. However, it is typically related to the growth in population over a single generation or the time taken for eggs to hatch into larvae which would be something between a couple of days and a few weeks. Parasitoids *P*_{1} and *P*_{2} parasitize uninfected hosts at intrinsic rates *α*_{1} and *α*_{2}, respectively. Handling and refractory times are modelled using a function of the form, (1−e^{−aH}), where *H* is the host population and *a* is a parameter which determines the efficiency of parasitoids and is related to the handling time associated with parasitizing hosts. This is known as an Ivlev functional response and is an alternative to the Holling type II functional response often used to model predation or host–parasitoid interactions (Sherratt *et al*. 1995; Pearce *et al*. 2006). Each parasitized host gives rise, on average, to *e*_{1} of the first parasitoid species *e*_{3} of the second in the next generation. Parasitoids *P*_{1} and *P*_{2} die at rates *d*_{1} and *d*_{2}, respectively. The mathematical model is, therefore,(2.1a)

(2.1b)

(2.1c)

There are four fixed points (steady states) in this system given byThe trivial fixed point (0, 0, 0) is always unstable. Any one of the other fixed points may be stable depending on the parameters. However, we note that an equilibrium point of coexistence does not exist.

In order to develop the model further, we now introduce host infection into the system. We denote uninfected and infected hosts *h*_{u} and *h*_{i}, respectively. We assume that contact between members of the host population is random and so we consider simple ‘mass action’ transmission with an infection term of the form *βh*_{i}*h*_{u} (DeJong *et al*. 1995; McCallum *et al*. 2001; Begon *et al*. 2002). As in the model without infection, hosts are subject to parasitism by two species of parasitoid *P*_{1} and *P*_{2}, the former laying a greater number of eggs in each clutch than the latter. Parasitism is once again modelled by the Ivlev functional response, (1–e^{–aH}), where *H*=*h*_{u}+*h*_{i} is the whole host population. Parasitoids *P*_{1} and *P*_{2} parasitize uninfected hosts at rates *α*_{1} and *α*_{3} and infected hosts at rates *α*_{2} and *α*_{4}, respectively. Each parasitized uninfected host gives rise, on average, to *e*_{1} of the first parasitoid species and to *e*_{3} of the second in the next generation. Parasitoids *P*_{1} and *P*_{2} die at rates *d*_{1} and *d*_{2}, respectively. We assume that the infection induces host mortality at a rate *γ*_{1}. If the host dies before the parasitoids emerge, then the juvenile parasitoids will also die. We assume that there is no recovery from the infection and that there is perfect vertical transmission. The model now including infection is, therefore,(2.2a)(2.2b)(2.2c)(2.2d)where *h*_{u} is the number of uninfected hosts, *h*_{i} is the number of infected hosts and *P*_{1} and *P*_{2} are the number of parasitoids of type 1 and 2, respectively, and *H*=*h*_{u}+*h*_{i}.

To non-dimensionalize the system, we set , , , and . Dropping the tildes gives the following non-dimensional system:(2.3a)(2.3b)(2.3c)(2.3d)where *s*_{i}=*α*_{i}/*r* *ρ*_{i}=*a*_{i}*K* *ν*=*Kβ*/*r* *c*_{i}=((*e*_{i}*α*_{i})/*r*) and *μ*_{i}=*d*_{i}/*r*.

We envisage a situation in which the parasitoids are gregarious, so the conversion rates *e*_{i} may be greater than 1 and if *e*_{i} is sufficiently high, it is quite possible for *c*_{i} to be greater than *s*_{i} in the non-dimensionalized system. If a host dies prior to parasitoid emergence, the parasitoids it is supporting will also die, so the conversion rate must be reduced to take into account host mortality. Given that infected hosts suffer increased mortality under stress, we require *e*_{1}>*e*_{2} and *e*_{3}>*e*_{4} and also that *c*_{1}>*c*_{2} and *c*_{3}>*c*_{4}. The parameters discussed below represent a situation where *P*_{2} has a lower conversion rate than *P*_{1} on uninfected hosts, but this conversion rate is less affected by the infection. Such a situation might arise through *P*_{2} laying smaller clutches of eggs and so that the host is subjected to lower stress levels. Disease can induce detectable physiological changes in hosts. It is not, therefore, unreasonable to suppose that a parasitoid might have different rates of parasitism of infected and uninfected hosts, because it may be easier to parasitize a weaker host, but the chance of the host dying before the juvenile parasitoids are ready to emerge is higher. Within this general framework, we consider three particular cases that illustrate the range of dynamics possible with the inclusion of infection in the system.

Case(i): using reference parameters *s*_{1}=0.85, *s*_{2}=0.6, *s*_{3}=0.7, *s*_{4}=0.88, *ρ*_{i}=2.5, *ν*=0.17, *γ*=0.16, *μ*_{i}=0.2, *c*_{1}=0.9, *c*_{2}=0.2, *c*_{3}=0.4 and *c*_{4}=0.35, we obtain nine biologically realistic fixed points:It should be noted that in the absence of one species of parasitoid, there is no fixed point of coexistence—either infected or uninfected hosts become extinct and this is the case regardless of parameter choice. Similarly, if the infection is omitted from the system, the weaker parasitoid competitor invariably goes extinct and, again, there is no coexistence fixed point. Both species of parasitoid and the infection are, therefore, necessary for coexistence despite the fact that the infection is detrimental to all individuals.

Linear stability analysis shows that with the exception of , which is a stable spiral (as illustrated in figure 1), the fixed points are all unstable. Numerical simulations show that initially the system exhibits high-amplitude, high-frequency oscillations which decrease in size as the host infection levels rise. After a peak in the level of infection, we see damped oscillations and the system rapidly reaches a steady state as the populations are attracted to the coexistent fixed point. By altering the initial conditions, it is possible to extend the length of time the transient dynamics persist for as can be seen from figure 1*b*.

Case(ii): we now vary the parameters *ν*, *μ*_{1} and *μ*_{2} and set the parameters to *s*_{1}=0.85, *s*_{2}=0.6, *s*_{3}=0.7, *s*_{4}=0.88, *ρ*_{i}=2.5, *ν*=0.15, *γ*=0.16, *μ*_{1}=0.16, *μ*_{2}=0.085, *c*_{1}=0.9, *c*_{2}=0.2, *c*_{3}=0.4 and *c*_{4}=0.35 as in figure 2. The fixed point of coexistence is now an unstable spiral and a limit cycle exists. Again the system initially exhibits the transient high-amplitude, high-frequency oscillations but in this case, when the infection levels rise, although the oscillations decrease in size, they do not slow down and soon the system converges to a limit cycle about the coexistent fixed point (see figure 2*c*).

Case(iii): if we now consider simply increasing the infection rate *ν* so that the parameters are *s*_{1}=0.85, *s*_{2}=0.6, *s*_{3}=0.7, *s*_{4}=0.88, *ρ*_{i}=2.5, *ν*=0.2, *γ*=0.16, *μ*_{i}=0.2, *c*_{1}=0.9, *c*_{2}=0.2, *c*_{3}=0.4 and *c*_{4}=0.35, then numerical simulations of the system produce the dynamics shown in figure 3*a*. The only difference between this set of parameter values and the one which gave the stable spiral in figure 1 is an increased rate of infection *ν* between hosts. However, the dynamics are very different. We see pulses of high-frequency, large amplitude oscillations which vary in length and regularity. Even after 50 000 time-steps, the system has not settled down (see figure 3*b*). Linear stability analysis reveals that the fixed point has associated eigenvalues −0.065±0.274i and 0.002±0.031i. Locally, therefore, this fixed point is unstable (a saddle point in ^{4}). However, all the other fixed points are also unstable, but the solutions remain bounded. A phase plot of *h*_{u}, *h*_{i} and *P*_{1} shown in figure 3*c* does not suggest that the solution trajectories are approaching any kind of a limit cycle, indeed the orbit looks chaotic. The normal test for chaos is to calculate the largest Lyapunov exponent *λ* of a system (Gottwald & Melbourne 2004). If *λ*>0, then nearby trajectories diverge exponentially, whereas if *λ*<0, then nearby trajectories stay close together. Thus, positive *λ* implies that a system is chaotic and negative *λ* that it is not. The Lyapunov exponent of the system as calculated by MATDS using the algorithm in (Wolf *et al*. 1985) is 0.05544>0. There is much discussion about what constitutes chaos and no computer simulation can show that a system is not periodic at a range just outside the length of the simulation. What we can say is that the system is not periodic within 50 000 time-steps—a time-scale more associated with evolution than population dynamics—and within that time we see chaos-like behaviour.

## 3. Spatio-temporal dynamics

Parasitoid–host systems in nature exist in a complex mosaic of spatial and temporal heterogeneity, and there has been an enormous volume of experimental work, both field and laboratory, which attempts to understand and describe the way in which parasitoids respond to this patchiness in their resource. In general, it is recognized that spatio-temporal heterogeneity can impart considerable stability to systems which might otherwise be unstable. The model introduced in §2 does not consider how spatial effects may influence the system. In reality, hosts and parasitoids are likely to move within a given spatial domain, with individuals encountering one another spatially. In order to model this more realistic setting, we therefore now allow each species to move randomly throughout a given domain by including random motility (diffusion) in the model. In some systems, parasitoids are known to use chemical cues (kairomones) in searching for hosts which enable them to aggregate at locations of high kairomone concentration (Vinson 1976; Schofield *et al*. 2002, 2005*a*). However, we do not consider such a situation here. For simplicity, we consider a one-dimensional domain *Ω*=(0, *L*) and impose zero-flux boundary conditions on ∂*Ω* to close the system. We assume that the uninfected and infected hosts have random motility coefficients *D*_{hu} and *D*_{hi}, respectively, while the species of parasitoid *P*_{1} and *P*_{2} have random motility coefficients *D*_{p1} and *D*_{p2}, respectively. We set and use the same substitutions as in obtaining system 3 to obtain the non-dimensional partial differential equation model below.(3.1a)(3.1b)(3.1c)(3.1d)where , , and and all other parameters are as before.

Throughout this section, we take *D*_{1}=*D*_{2}=8×10^{−7} and *D*_{3}=*D*_{4}=7.5×10^{−6}. We have set host diffusion lower than parasitoid diffusion because transmission of infection generally occurs between the larval stages and wasps move faster than host larvae. However, we still observe the same dynamics when all diffusion coefficients are equal. The initial conditions are and (figure 4*a*).

Figure 4 shows snapshots of the (numerical) solutions of the PDE system with underlying ODE dynamics given in figure 1 where the fixed point is a stable spiral. Initially, a wave of uninfected hosts closely followed waves of *P*_{1}, and at lower levels *P*_{2} leads a wave train through the domain (figure 4*b*). However, in the absence of infection *P*_{2} cannot compete and drops to low levels. We then see rich spatio-temporal dynamics associated with the transient phase from figure 1 with high-frequency, large amplitude oscillations (figure 4*c*). As the host infection invades, *P*_{2} follows it and the oscillations slow down and rapidly decrease in amplitude (figure 4*d*,*e*) and by *t*=2500, the system has reached its steady state (figure 4*f*). The transient oscillations in the healthy host population mean that the infection is not invading a homogeneous population and the wave of invasion is disrupted. The order of invasion can be explained by the fact that parasitoids require hosts to reproduce and hence cannot exist without them, so the hosts must invade first. The infection invades at a slower speed than both the healthy hosts *h*_{u} and the parasitoids *P*_{1}. Similarly, *P*_{2} has a very low rate of population growth in the absence of infection and it does not invade until the infection has established itself.

Figure 5 shows snapshots of the (numerical) solutions of our PDE system corresponding to the underlying ODE dynamics of figure 2. Here again, we see a travelling wave of healthy hosts followed closely by a travelling wave of *P*_{1} and *P*_{2} (figure 5*a*). Again, *P*_{2} cannot compete and we see a much longer period of transient dynamics (figure 5*b*). The disruption of the wave of infection is extreme here. It does not follow a wave-front, but reaches low levels throughout the domain before rising to a peak (figure 5*c*) and spreading out from it. *P*_{2} cannot compete in the absence of infection, so its population density increases behind the peak of infection and it takes much longer to become established (not until *t*=1800 in figure 5*d*). The oscillations decrease in amplitude (figure 5*e*), but this is temporary and the amplitude increases slightly again. The long-term dynamics of the underlying ODE system are a stable limit cycle, and correspondingly in the PDE solutions, we see persistent rich spatio-temporal heterogeneity as shown in figure 5*f*. This type of dynamic spatio-temporal behaviour has been observed in previous work examining predator–prey systems and host–parasitoid systems (cf. Sherratt *et al*. 1995; Pearce *et al*. 2006).

Figure 6 shows snapshots of the (numerical) solutions of our PDE system corresponding to the underlying (chaotic) ODE dynamics of figure 3. The speed of the wave of infection is faster here because the rate of infection is higher. As a result, the population of parasitoid species 2, *P*_{2}, becomes established throughout the domain much more quickly. The peaks in host infection seen in figure 3*a*,*b* do not occur when diffusion is introduced to the system, but we do see pulses of fast oscillations and as can be seen from figure 6*f*, the infection level continues to vary widely both spatially and temporally. The complex spatio-temporal dynamics seen in figure 5 are very much in evidence here with larger oscillations and a greater degree of spatio-temporal heterogeneity.

## 4. Discussion

The mathematical model presented in this paper develops a theoretical framework with which to address the issue of infection in host–multi-parasitoid systems. It combines two well understood mechanisms—simple contact transmission of infection and host-mediated parasitoid competition. The results of our model have shown that if one species of parasitoid is removed from the system, there are no coexistent fixed points and that the inclusion of multiple parasitoids allows for persistence of the disease. On the other hand, without the infection, competitive exclusion means that either one or both species of parasitoid will inevitably become extinct. The addition of disease dynamics to the model, therefore, allows for coexistence of all species, so although the disease is detrimental to all agents in the system, it indirectly benefits the less efficient competitor. The prediction that the two parasitoid species can only exist in the presence of the pathogen is an interesting and novel result, because it represents a departure from the predictions of some models of apparent competition in host–parasitoid systems, where a shared host can be detrimental to the survival of one of the parasitoid species (Bonsall *et al*. 2002). The operation of this process can be seen in figure 2, where pulses of infection reduce the population of the first parasitoid species, which are then followed by smaller peaks in the population of the second parasitoid. It also differs from classic two-consumer–two-resource models such as those discussed in Vandermeer & Pascual (2006), because the interactions between the healthy and infected host populations are not simply based upon resource competition. If the parasitoids were removed from the system, we would see a classic susceptible–infected disease model.

Our model generates chaos-like dynamics for a wide range of parameter values, but even when this does not occur, the system initially exhibits complex transient dynamics which can be extended by simple alteration of the initial conditions. These transient dynamics remain important when one-dimensional diffusion is included in the model. The existence of these transient chaos-like dynamics suggests that in any system which is subject to regular substantial disturbance, especially where the initial starting conditions are repeatedly re-established, dynamics which are essentially chaotic could persist almost indefinitely. Sait *et al*. (2000) observed that transient dynamics persisted in an experimental host–parasitoid–pathogen system for ecologically significant periods, but their experiments used only one parasitoid species, and invariably the system became extinct. The parameters for the three case studies in this paper have been chosen to illustrate the rich variety of dynamics possible from a simple model. In all cases in this paper, we note that *c*_{1}>*c*_{3} and *c*_{4}>*c*_{2}. However, if the *ρ*_{i} are allowed to differ from each other, then coexistence can still occur with *c*_{1}>*c*_{3} and *c*_{2}>*c*_{4}. Finally, we note that in this paper we have assumed that vertical transmission of the infection is perfect and that parasitoid movement is a random process. Further development of this work will include the effect of varying vertical transmission rates and directed parasitoid movement. Many parasitoids use both chemical and visual cues when searching for hosts (Vinson 1976; Schofield *et al*. 2005*a*,*b*) and a model which incorporates aggregation behaviour via chemotaxis is currently under investigation.

Agro-ecosystems are examples of semi-natural habitats where there is regular, often seasonally based, disturbance which returns them to their initial state or something very close to it. If a system were subject to such continual small perturbations, it could exhibit ongoing transient dynamics and never arrive at a steady state or periodic limit cycle, even though it had the innate capacity to do so. This suggests that prescriptions designed to promote biodiversity in repeatedly disturbed systems, which are common in UK agriculture and the natural areas surrounding them, for example, may be less successful than hoped simply owing to the regular disturbance. This might be especially detrimental in the case of agricultural systems, because the perturbation is always negative, in the sense that all populations within the association are likely to be reduced considerably at each intervention. In addition, any benefits which such prescriptions might have in enhancing the effects of natural enemies on pest populations within the crop by acting as reservoirs is likely to be similarly reduced. Although we have discussed host–parasitoid interactions, the model could equally be applied to predator–prey–pathogen interactions (Hassell 1978) and, more generally, our model suggests that in any situation where multiple infections occur, removal of one detrimental (to the host) agent may lead to the loss of other detrimental agents in the association. This may have implications for the study of disease dynamics and the treatment of pathogenic infections.

## Footnotes

- Received September 26, 2006.
- Accepted November 1, 2006.

- © 2006 The Royal Society