William E. Snyder and Anthony R. Ives. 2003. Interactions between specialist and generalist natural enemies: parasitoids, predators, and pea aphid biocontrol. Ecology 84:91–107.

We constructed a model to investigate the interactions between parasitoids and predators, and how these interactions affect pea aphid control. We have four interrelated objectives. First, field experiment 3 demonstrated additivity between the effects of parasitoid and predators. Nonetheless, mortality from the fungal epizootic was high. In the absence of the fungus, would parasitism and predation be additive? Second,  field experiment 3 was run for 21 d, at which time the aphid densities within cages had peaked.  The duration of the experiment was therefore appropriate for determining the effect of biological control on peak aphid densities within a harvesting cycle. Nonetheless, the average interval between harvests is longer (roughly 32 d), and the population size of aphids at the end of the harvesting cycle will influence both the dispersal of aphids to other fields and the number of aphids at the start of the following harvesting cycle. How did the duration of experiment 3 affect our conclusion about the additivity of parasitism and predation? Third, although containing hundreds of stems, the cages were small compared to an entire field, and the cages excluded movement of predators in response to aphid densities. Would aphid density-dependent migration of predators affect the additivity between parasitism and predation? (We do not investigate the potential effects of A. ervi migration, because the dynamical coupling between aphids and parasitoids makes this problem much more complex [e.g., Ives 1992, Murdoch et al. 1992].) Fourth, although our experiments involve the specific system of pea aphids in alfalfa, we are interested in the broader question of whether additive interactions should be expected in other host–parasitoid–predator systems. This requires understanding the underlying processes leading to additivity. What potential processes could lead to nonadditivity in our aphid-parasitoid-predator system and other similar systems?

Model structure. – Because we want to fit the model to our experimental data, we have tried to strike a balance between including enough details to describe aphid and natural enemy dynamics and excluding enough details to limit the number of parameters that need to be estimated. The model explicitly includes stage structure of the aphid population, dividing the population into the five aphid instars. Thus, the model incorporates the development time of aphids and also variation among aphid instars in their susceptibility to parasitism. To account for the development time of parasitoid larvae within aphids, the parasitoid population is also treated as stage structured, with stages corresponding to the developmental time of aphid instars. We do not model predator dynamics explicitly, but instead implicitly incorporate predation as a decline in the stage-specific survival of aphids (and larval parasitoids within still-living aphids). Because we have no information about the mechanisms of transmission and infection of the fungal pathogen from our experiments, we treat the epizootic as a delayed density-dependent source of mortality on aphids. Delayed density-dependent aphid population growth could also be caused by aphid density-dependent declines in plant quality, and we make no attempt to separate this from the effects of the fungal pathogen on aphid dynamics. We fit the model to the data on aphid and mummy densities from all 3 field experiments separately; we did not include the data on fungal cadavers in model fitting, because the pathogen is incorporated implicitly rather than explicitly in the model.

The aphid dynamics are governed by three processes that are iterated on a 2-d time scale: density-independent population growth, density-dependent survival, and parasitism. The 2-d period of iteration corresponds roughly to the development time in the four juvenile instars, creating a 10-d period from first instar to reproductive adult. The density-independent population growth is given by the Leslie matrix, L, with the form

(A.1)

where s_a is the survival of aphids during an instar and f is the fecundity of adults in a 2-d period. For simplicity, we assume that survival in all instars is the same, recognizing that in reality this is not the case. Letting X(t) denote the vector of densities of aphids in each of 5 instars, the densities following reproduction and density-independent survival, X'(t), are given by

We assume that density-dependent survival of aphids occurs with a time delay. A time delay is appropriate to model both the density-dependent response of the fungal pathogen to increases in aphid density and a density-dependent decline in plant quality that occurs as plants succumb to high aphid loads. In fitting the model to the experimental data, we investigated density-dependent time delays of 0, 2, 4, and 6 d. The best-fitting model with a 4-d time delay fit the data better than the best-fitting models with 0-, 2-, and 6-d time delays, so we used a 4-d time delay. Thus, if X''(t) denotes the vector of aphid instar densities following density-dependent survival,

where x_s(t - 2) is the total aphid density (sum of all instars, including parasitized aphids) 2 time steps (4 d) preceding t, and k is a parameter governing the strength of density dependence; when k = 0, there is no density dependence, while increasing k corresponds to increasing density dependence. We assume that all instars experience the same strength of density-dependent survival.

Parasitoids are assumed to have a functional response governed by a parameter g; when g = 0, the functional response is Type I, while larger values of g correspond to greater parasitoid saturation leading to stronger Type II functional responses. Observational studies on the functional response of A. ervi (Ives et al. 1999) have demonstrated that this parasitoid has a functional response close to Type I.  Nonetheless, we used a functional form permitting a Type II functional response in the model both to confirm the Type I functional response of A. ervi and to allow us to consider the consequences of a Type II functional response in subsequent analyses. Furthermore, the model incorporates instar-specific differences in susceptibility to parasitism as revealed through behavioral studies (Ives et al. 1999). If x_i''(t) and x_i(t + 1) denote the densities of unparasitized aphids in instar i before and after parasitism, respectively, and y₈(t) is the density of parasitoid adults, then

where a is the overall searching efficiency of parasitoids and r_i is the relative searching efficiency on instar i.  For values of r_i, we use those obtained from studies of A. ervi: 0.12, 0.27, 0.39, 0.16, and 0.06 for instars 1 – 5, respectively (Ives et al. 1999).

To account for parasitoid developmental time, the parasitoid population is divided into 8 stages. The first 4 stages (8 d) are for larvae within living aphids, the following 3 stages (6 d) are for mummies, and the final stage is adults. The number of parasitoids entering the first stage is determined by the parasitism rate on aphids. Thus, if y_i(t + 1) denotes the density of parasitoids in stage 1 at time step t + 1,

(A.5)

Because parasitoids exist within living aphids, they experience the mortality that affects aphids. Following the experimental results of Rauwald and Ives (2001), we set the survival of parasitized aphids equal to that of unparasitized aphids. Thus, for stages i = 1 – 4,

We assume that mummies have density-independent survival s_m, so for stages i = 5 – 6

Finally, we assume that adult parasitoids have survival s_w and do not senesce, so

This model implicitly includes generalist predators in the parameters s_a and s_m. Increasing the density of predators will decrease the density-independent survival of aphids (s_a) and/or mummies (s_m). The model implicitly includes the fungal pathogen and possible delayed density-dependent effects of aphid density on plant quality through the parameter k.

Parameter estimation. – We fit the model separately to each of the three field experiments. We estimated parameters under the assumption that all differences between model and data were due to measurement error, thereby allowing a total least-squares fitting procedure (Ludwig and Walters 1989, Hilborn and Walters 1992). For each cage in a given experiment, we generated simulated data by starting with initial estimates of aphid and parasitoid densities and iterating the model for the equivalent length of the experiment. The fit of the simulated data to the real data was assessed by the sum of squared differences between the log aphid densities and log mummy densities. Comparing differences in log densities rather than absolute densities makes the assumption that the coefficient of variation of the measurement error is independent of mean aphid and parasitoid densities. Furthermore, since aphids and mummies were sampled in the same way, we assumed that the coefficient of variation of their measurement errors were the same. Using the sum of squared differences in log densities to quantify the model fit, we estimated parameters by selecting those values that gave the lowest sum of squared differences using a numerical minimization routine (the function fmins.m in Matlab, Mathworks 1996).

The fit of the model to data proved to be very sensitive to the initial values of aphid density and percent parasitism. In experiment 1, we treated initial aphid density in each cage as parameters to be estimated by the model (11 total parameters, one per cage). We used two parameters for the initial level of parasitism of aphids, one parameter for the treatments that received supplemental aphids and one for treatments with ambient initial aphid densities; this accounts for the fact that the aphids used to supplement cages were unparasitized. Finally, we scaled the densities of the introduced adult parasitoids using two parameters to separate the cages with and without supplemental aphids; preliminary analyses showed that the relative impact of introduced parasitoids was greater in the cages with ambient aphid densities than in cages with supplemented aphids. Initial mummy densities were not estimated, but instead observed initial values were used. In experiment 2, for starting aphid densities we used the observed values in all cages scaled by a single parameter. This scaling parameter was used because the aphids in the first census were sampled by D-vaccing an entire cage, while subsequent samples were made by counting aphids on stems. Preliminary analyses revealed that aphid densities estimated from D-vaccing were roughly double those estimated from stem counts. For experiment 3, we used the observed aphid densities as the initial conditions. For initial parasitism, we used a single parameter that was the same for all cages.

Because the model includes stage structure of both aphids and parasitoids, initial stage structures must be set for the simulations. For all experiments, we assumed that the initial stage structure of unparasitized aphids was given by the stable age distribution computed by the density-independent Leslie matrix (Eq. A.1). The stage distribution of parasitized aphids was set equal to that of the unparasitized aphids, while parasitoid larvae were assigned equally among parasitoid stages. The delayed density-dependence in aphid survival (Eq. A.3) requires values of total aphid density for 2 and 4 d preceding the initial sample. For these values, we use the Leslie matrix (Eq. A.1) to backcast the density of aphids that would occur if the population was growing exponentially at its stable-stage distribution.

Two sets of parameters were difficult to fit simultaneously. First, estimates of aphid fecundity, f, and aphid survival, s_a, strongly covaried, with the statistical fitting procedure poorly distinguishing between the case of high fecundity and low survival vs. low fecundity and high survival. This is a common problem in fitting stage-structured models (Wood 1994). To overcome this problem, for all experiments we fixed fecundity f = 8, which is the maximum 2-d fecundity for pea aphids obtained in laboratory experiments (Thiboldeaux 1986). Second, estimates of adult parasitoid survival, parasitoid searching efficiency, and mummy survival all strongly covaried; the statistical fitting procedure poorly distinguished between the case of low adult survival, high searching efficiency and low mummy survival vs. high adult survival, low searching efficiency and high mummy survival. To solve this, we set adult parasitoid survival, s_w, to 0.8, and then fit the other two parameters. Because experiment 1 had much greater differences among cages in aphid and parasitoid dynamics than experiment 3, we estimated the searching efficiency and mummy survival from experiment 1, and used the resulting estimate of searching efficiency in experiment 3 (experiment 2 had no parasitoids).

Our intention in fitting the model to our field data was to give a tool to explore how interactions among natural enemies influence pea aphid control. This requires a model that can mimic the population dynamics observed in the experiments. Because the model makes numerous simplifying assumptions, we have less confidence in the values of the estimated parameters. For example, the model assumes that all aphid stages last 2 d and that stage-specific aphid survival is the same for all stages. Because these assumptions are only approximately correct, we do not have great faith in the precise estimate of aphid survival, s_a. Due to the uncertainty in parameter estimates caused by uncertainty about the structure of the model, we do not attempt to obtain confidence limits for the estimates of the parameters.

Literature cited

Hilborn, R. and C. J. Walters. 1992. Quantitative fisheries stock assessment. Chapman and Hall, New York, New York, USA.

Ives, A. R. 1992.  Continuous-time models of host–parasitoid interactions.  American Naturalist 140:1–29.

Ives, A. R., S. S. Schooler, V. J. Jagar, S. E. Knuteson, M. Grbic, and W. H. Settle. 1999.  Variability and parasitoid foraging efficiency: a case study of pea aphids and Aphidius ervi.  American Naturalist 154:652–673.

Ludwig, D., and C. J. Walters. 1989. A robust method for parameter estimation from catch and effort data. Canadian Journal of Fisheries and Aquatic Sciences 46:137–144.

MathWorks, I. 1996. MATLAB. The MathWorks, Inc., Natick, Massachussetts, USA.

Murdoch, W. W., C. J. Briggs, R. M. Nisbet, and W. S. C. Gurney. 1992.  Aggregation and stability in metapopulation models.  American Naturalist 140:41–58.

Rauwald, K. S., and A. R. Ives. 2001.  Biological control in disturbed agricultural systems and the rapid re-establishment of parasitoids.  Ecological Applications 11:1224–1234.

Thiboldeaux, R. 1986. The effect of temperature on population level interactions between the pea aphid, Acyrthosiphon pisum, and its primary parasitoid, Aphidius ervi.  Thesis. University of Wisconsin, Madison, Wisconsin, USA.

Wood, S. N. 1994. Obtaining birth and mortality patterns from structured population trajectories. Ecological Monographs 64:23–44.