Comparing predator-prey models to Luckinbill’s experiment with Didinium and Paramecium
Gary W. Harrison
Editor’s Note
Progress in science is supposed to depend upon interaction between theory and experiment, but often this interaction is too weak to be effective. Harrison’s paper illustrates the gains that are possible when both the data and the theory are given serious attention. Harrison confronts a variety of predator-prey models with Luckinbill’s (1973) data on Didinium and Paramecium. A standard model gives qualitative agreement: this is the point where too many of us declare a victory and turn to other pursuits. Harrison is not satisfied, but goes on to fit a variety of models to the data and examine the deviations between the data and the fits in order to understand the models and suggest improvements. Few models can stand up to such treatment. The apparent winner is a new model that allows for a lag between consumption and reproduction. This model also provides good fits when the experimental medium is thickened. This investigation may provide insight for further experiments on the details of the predator-prey interaction.
Introduction
Mathematical models of predator-prey interactions, starting with the models of Lotka (1925) and Volterra (1931) and including the many refinements and modifications of their model (e.g., Murdock and Oaten 1975, Berryman 1992), form a basis for the theory of predator-prey interactions, but have not been tested often enough against actual experimental data. The purpose of this paper is to test standard predator-prey models against the data from an experiment of Leo Luckinbill (1973) to see how well the predictions of the models match the data both qualitatively and quantitatively and to determine what refinements or modifications are most important to obtain a good match.
All models show that theoretical predator-prey systems have a tendency to oscillate, but researchers have had a difficult time duplicating these results in laboratory settings because usually the predators quickly consume all the prey and then starve. Luckinbill (1973) grew Paramecium aurelia together with its predator Didinium nasutum, in 6 mL of cerophyl, which provides a food source for the Paramecium. The Didinium consumed all of the Paramecium within a few hours [ILLUSTRATION FOR FIGURE 1 OMITTED], similar to the results of Gause (1934), but when the solution was thickened by adding methyl cellulose, both populations went through two or three cycles over a period of 10-18 d before the Didinium became extinct [ILLUSTRATION FOR FIGURE 2 OMITTED]. The methyl cellulose slowed the movement of both species, but did not add any physical refuge for the prey as in Gause et al. (1936) and Flanders (1948) nor any physical complexity as in Huffaker (1958); both species were free to move throughout the medium. When the medium was prepared with methyl cellulose and half-strength cerophyl, the populations maintained sustained oscillations for 33 d, at which time the experiment was terminated [ILLUSTRATION FOR FIGURE 3 OMITTED]. How well can current predator-prey models explain these results of Luckinbill’s experiment?
It is now generally accepted that a realistic model of predator-prey interaction should include a carrying-capacity for the prey and a satiation effect in the predator functional response. The section Qualitative description below shows that adding methyl cellulose should reduce the half-saturation constant in the functional response and using half-strength cerophyl should reduce the carrying capacity, and these two features are sufficient to qualitatively explain Luckinbill’s results. The quantitative fit between this basic model and the observed population densities is not completely satisfying, however. Therefore, various refinements – such as predator intraspecific competition, Leslie type models, “ratio-dependent” predation rates, sigmoid functional response, and time delays in the response of reproductive rates to increased food – are tested to see which gives the best fit to the observed data. Four variations in the error criteria were used, but produced only minor differences in the goodness of fit. Predator mutual interference – combined with a predator reproductive rate that depends on nutrient and/or energy storage in the predator population and thus, in a delayed manner, on the consumption rate – gives an excellent fit, but the best model is a combination of the delayed predator reproductive rate with a sigmoid functional response.
The models make two predictions that could be tested experimentally. First, further reduction of the food supply for the Paramecium and/or thickening of the medium could produce a nonoscillating stable steady state for the population densities. Second, the best model requires’ a more sigmoid functional response to fit the data in the thickened than in the unthickened medium, indicating that not only is the overall effective rate of search by the Didinium lower in the thickened medium, but also there is a more pronounced reduction in their search rate at low prey densities when the effort of swimming is greater.
Qualitative Description
Basic model
Let x be the prey density (numbers per cubic centimetre) and y the predator density. The basic predator-prey model is
dx/dt = p(1 – x/K)x – [Omega]f(x)y (1) dy/dt = [Sigma]f(x)y – [Gamma]y. (2)
The parameter [Rho] is the specific growth rate of the prey, K is the carrying capacity of the prey in the absence of predators, and [Omega]f(x)y is the predation rate when the prey have density x and the predators have density y. The rate of prey consumed per predator, [Omega]f(x), is called the functional response of the predator. There is good evidence (Holling 1959, Murdock and Oaten 1975) that many predators exhibit a saturation effect in their functional response when prey is abundant, that is
[Omega]f(x) = [Omega]x/[Phi] + x, (3)
which is known as a Holling Type 2 functional response. Here [Omega] is the maximum rate of prey consumption by a single predator and [Phi] is the half-saturation constant, the level of prey at which half the maximum consumption rate occurs. Eq. 2 follows the usual assumption that the predator rate of increase (numerical response) is proportional to the rate of prey consumption, with the efficiency of converting prey to predator given by [Sigma]/[Omega]. There is less experimental evidence for the proper form of the numerical response, but Slobodkin (1986), Beddington et al. (1976), and Coe et al. (1976) give evidence that the proportionality assumption used here is valid. The term – [Gamma]y represents the rate at which the predators would decrease in the absence of the prey.
Although differential equation predator-prey models such as Eqs. 1-3 are sometimes considered too simplistic to give realistic results, this model can give a qualitative explanation of Luckinbill’s experimental results. To understand how, one must first review the stability theory of this model.
Qualitative stability theory. – The stability theory of the basic model in Eqs. 1 and 2 with functional response (Eq. 3) is well understood. The prey isoclines where dx/dt – 0 are x = 0 and
y = [Rho]/[Omega](1 – x/K)([Phi] + x), (4a)
which is a concave-down parabola with vertex at x = (K – [Phi])/2, and x intercept at x = K. The predator isoclines where dy/dt = 0 are y = 0 and the vertical line
x = [Gamma][Phi]/[Sigma] – [Gamma]. (4b)
Thus if [Gamma][Phi]/([Sigma] – [Gamma]) [less than] K there is an equilibrium
[E.sub.1] = ([x.sup.*], [y.sup.*]),
with
[x.sup.*] = [Gamma][Phi]/[Sigma] – [Gamma], [y.sup.*] = [Rho]/[Omega](1 – [x.sup.*]/K)([Phi] + [x.sup.*]), (4c)
that has both species present. The graphical stability condition of Rosenzwieg and MacArthur (1963), which was proved algebraically by Freedman (1976), is that [E.sub.1] is stable if the vertical isocline is to the right of the vertex of the prey isocline, that is if
[Gamma][Phi]/[Sigma] – [Gamma] [greater than] K – [Phi]/2. (5)
(Throughout this paper, stable will mean asymptotically stable in the strict mathematical definition.) Harrison (1986) gives a more complete exposition of these results. When [E.sub.1] is stable, any oscillations in predator and prey densities will eventually be damped out and both species will coexist at the levels [x.sup.*] and [y.sup.*].
If Eq. 5 is not satisfied, then [E.sub.1] is unstable, and a trajectory that starts close to [E.sub.1] will exhibit diverging oscillations that spiral out toward a stable limit cycle, whose existence can be shown by using Poincare Bendixson theory (Hastings 1978, Hsu 1978). Numerical solutions show that the further [E.sub.1] is to the left of the vertex of the prey isocline, the faster the oscillations diverge and the larger the amplitude of the limit cycle. Once the trajectory is close to the limit cycle, the model will exhibit continued oscillations, but whether the real predator and prey can coexist depends on how close these oscillations come to the boundaries x = 0 and y = 0. A density of [less than]1/6 individuals/mL in 6 mL of liquid implies extinction in the real system, even though the mathematical model treats any density greater than zero as a viable population. Note, however, that in 600 mL of liquid a density of 1/6 individuals/mL implies 100 individuals alive, which is one reason that natural populations can survive even though small experimental populations go extinct.
The effect of adding methyl cellulose. – The relation between the stability theory and the results of Luckinbill’s experiment depends on how the addition of the methyl cellulose would change the parameters in Eqs. 1-3, particularly the coefficients [Phi] and [Omega] in Eq. 3. The interpretation of [Phi] and [Omega] is clarified by examining Holling’s (1959) formula for the Type 2 functional response:
[Mathematical Expression Omitted],
where h is the handling time for each prey captured and e is the rate of effective search by the predators. Thus the parameters in Eq. 3 for the maximum predation rate [Omega] and the half-saturation constant [Phi] are given by
[Omega] = 1/h; [Phi] = 1/eh. (6)
The primary effect of adding methyl cellulose to the experimental medium was to reduce the swimming speed of both species, which would reduce the searching effectiveness e of the Didinium, and according to Eq. 8 would increase [Phi]. Any effect on the handling time and hence on [Omega] would be minimal, as would the effect on any of the other parameters in Eqs. 1 and 2.
In the original mixture the searching effectiveness of the Didinium was great enough that [Phi] must be very small, with the result that [x.sup.*] given by Eq. 4c is very small and the equilibrium is far to the left of the vertex in the prey isocline. From Inequality 5 the equilibrium is unstable and the diverging oscillations around it quickly carry densities to such low levels that no Paramecium are left [ILLUSTRATION FOR FIGURE 4a OMITTED]. Adding the methyl cellulose dramatically reduces the searching effectiveness of the Didinium so that [Phi] must be much larger. This makes [x.sup.*] much larger, but apparently still not large enough to be to the right of the vertex in the prey isocline. The equilibrium is still unstable, but the diverging oscillations are centered about a point further from the axis and diverge more slowly, so that the system goes through three cycles before the density of either species becomes low enough to represent extinction in 6 mL of liquid [ILLUSTRATION FOR FIGURE 4b OMITTED].
The effect of using half-strength cerophyl mixture. – Reducing the nutrient supply for the Paramecium may reduce the specific growth rate [Rho] somewhat but would most dramatically reduce the carrying capacity K. Indeed, Luckinbill’s experimental results when the Paramecium was grown by itself show that the carrying capacity was reduced from [approximately equal to]900 individuals/mL to [approximately equal to]400 individuals/mL. Reducing K has no effect on the vertical predator isocline, but shifts the prey isocline (Eq. 4a) down and to the left [ILLUSTRATION FOR FIGURE 4c OMITTED]. Thus it is much more likely that Inequality 5 will be satisfied and the equilibrium point [E.sub.1] is stable, or even if Inequality 5 is not satisfied, [E.sub.1] is much closer to the vertex of the prey isocline and the stable limit cycle surrounding it will have smaller amplitude and be further from the axes. This is Rosenswieg’s “paradox of enrichment” (Rosenswieg 1969) in reverse. Reducing the nutrient supply of the prey is stabilizing, although in Luckinbill’s experiment it was not reduced enough to produce a stable equilibrium with constant population [TABULAR DATA FOR TABLE 1 OMITTED] densities but only a stable limit cycle with a small enough amplitude to avoid extinction.
Thus on a qualitative level the model can explain why both thickening the medium and reducing the nutrients had a stabilizing effect on Luckinbill’s experimental system. The model predicts that using even more methyl cellulose and/or even weaker cerophyl medium could produce a stable system with constant population densities [ILLUSTRATION FOR FIGURE 4d OMITTED].
Quantitative Fit of the Experimental Data: Methods
Although the model given by Eqs. 1-3 explains the experimental results on a qualitative level, it cannot give a very good quantitative fit to the experimental population densities reported by Luckinbill [ILLUSTRATION FOR FIGURE 5 OMITTED]. Therefore various modifications of the model were explored to see which could give the best fit, thus testing what other features may be important to include in models of predator-prey interactions. The results are summarized in Table 1.
Methods. – The Marquardt-Levenberg method was used to estimate the parameters in the differential equations that would give the least squares best fit to the experimental data, that is to find the parameters that would minimize [S.sup.2], defined by
[S.sup.2] = [summation of] [[x([t.sub.i]) – [x.sub.observed]([t.sub.i])/[s.sub.xi]].sup.2]
+ [[y([t.sub.i]) – [y.sub.observed]([t.sub.i])/[s.sub.yi]].sup.2], (7)
where x([t.sub.i]) and y([t.sub.i]) are the values computed from the differential equation, the observation times [t.sub.i] are at half-day intervals, N is the total number of observations, and [s.sub.xi] and [s.sub.yi] are weights representing the relative measurement errors in the observations, which are discussed below. The Marquardt-Levenberg method was implemented using the program given in Press et al. (1986) with some slight modifications. The differential equations were solved using a variable step size Runge-Kutta method, Felberg’s RK4(5) method (Danby 1985). The values of the empirical populations densities, [x.sub.observed]([t.sub.i]) and [y.sub.observed]([t.sub.i]), were taken from the graphs published by Luckinbill for Figs. 1 and 2 using a digitizer, which introduces another source of experimental error. For the final experiment [ILLUSTRATION FOR FIGURE 3 OMITTED], however, the data was provided directly by Leo Luckinbill (personal communication).
The data of Fig. 1 (without methyl cellulose) were not sufficient to effectively estimate all the parameters of the model. There were consistent estimates for [Omega]/[Phi] but not for [Omega] and [Phi] separately, for example. The data of Fig. 3 (half-strength nutrients) also proved very difficult to fit. A trajectory that looks like a close fit visually with the right number and amplitude of oscillations but which is slightly out of phase may produce a very large value for [S.sup.2]. Therefore the strategy adopted was to use the numerical method to find the best fit to Fig. 2, then as an added check to test whether varying [Phi] or K could produce reasonable fits to Figs. 1 or 3, respectively, without making any changes in the other parameters.
Error criteria. – Different values for the sum of squares error [S.sup.2] can be obtained depending on the weights [s.sub.xi] and [s.sub.yi] used in Eq. 7, the initial conditions used when solving the differential equations, and which parameters are estimated or held fixed. Four of these are reported in Table 1. The values [[S.sub.1].sup.2], [[S.sub.2].sup.2], and [[S.sub.3].sup.2] are all unweighted sum of squares errors, produced by letting [s.sub.xi] = [s.sub.yi] = 1 in Eq. 7. The computation of [[S.sub.2].sup.2], which keeps [Rho] fixed at 3.3, and [[S.sub.3].sup.2], which uses [t.sub.0] = 0.5 and estimates the initial conditions, are discussed below.
It may be argued that unweighted sum of squares error places too much emphasis on fitting the peaks at the expense of not fitting the sections of low densities. The other extreme is to assume that the observation errors are proportional to the densities observed, and measure relative errors by letting [s.sub.xi] = [x.sub.observed]([t.sub.i]) and [y.sub.observed]([t.sub.i]). This was tried but did not give good results. Since the errors at the peak densities were weighted very little, the algorithm selected parameters that produced flat, nonoscillating trajectories. A compromise between these two extremes is to let
[s.sub.xi] = 1 + [Alpha][x.sub.observed]([t.sub.i]); [s.sub.yi] = 1 [Alpha][y.sub.observed]([t.sub.i]). (8)
Experimenting with different values of [Alpha] showed that [Alpha] = 0.01 was a good choice, weighting the errors at the highest density about 1/6 as much as at the lowest density. The results are reported as [[S.sub.R].sup.2] in Table 1, but it produced only minor differences in the relative rankings of the different models. Fig. 5c shows that the fits at the low densities were not improved by using this relative error. Comparison of Fig. 5a, c with Fig. 9a, c shows that although the weighting makes a large difference in the best fit for the standard model, both weighting schemes produce best fit trajectories that appear almost identical for the best model found.
Table 1 also reports [R.sup.2], defined by
[R.sup.2] = SST – [[S.sub.1].sup.2]/SST,
[Mathematical Expression Omitted],
where [Mathematical Expression Omitted] and [Mathematical Expression Omitted] are the average values of the observed population densities and SST = total sum of squares. Since the differential equations are nonlinear and there is no evidence that the residuals are normally distributed, [R.sup.2] cannot be used to perform statistical tests, but it does give some indication of the fraction of the total variation explained by the model.
Reasonable parameter values. – It is necessary to make some a priori estimates of the parameters, both to give a starting point for the numerical procedure and to check whether the estimated parameters are reasonable. The plot of the Paramecium growing by itself [ILLUSTRATION FOR FIGURE 2 OMITTED] can be fit quite well by using a logistic equation (Eq. 1 with y = 0), with [Rho] = 3.3 (per day) and K = 898 (individuals). This value for K was used in all the models and was not one of the estimated parameters. Similarly, the Paramecium growing by themselves with reduced nutrients show a carrying capacity of [approximately equal to]400 [ILLUSTRATION FOR FIGURE 3 OMITTED], which was therefore used whenever a model was compared to the long data set shown in Fig. 3. The column [[S.sub.2].sup.2] in Table 1 reports the unweighted sum of squares error when p is held fixed at [Rho] = 3.3. In all the other cases [Rho] was estimated along with the other parameters, in which case a value of [Rho] that is not somewhere around [Rho] = 3 would cast doubt on the validity of the model. Comparing [[S.sub.1].sup.2] with [[S.sub.2].sup.2] shows that holding [Rho] = 3.3 makes [S.sup.2] slightly larger, but did not change any of the relative rankings of the models, although sometimes it had a large influence on the other parameter values.
A good estimate for [Omega], the maximum number of Paramecium captured per Didinium per day, is harder to make, but if one recalls that it is the inverse of the handling time (Eq. 6), it would not be unreasonable for [Omega] to be as big as 24 (one per hour) but probably not as big as 144 (one every 10 min). Wichterman (1986) quotes Reukauf (1930) that one Didinium can eat up to 12 Paramecium per day, which would imply [Omega] = 12, but Salt (1974) reports consumption rates up to two per hour, i.e., [Omega] = 48 Paramecium/d. One would expect [Sigma] in Eq. 2 to be smaller than [Omega] because of inefficiencies in converting prey into predator, but an individual Paramecium is larger than an individual Didinium (Wichterman 1986), so that [Sigma] may not be as small as expected. No attempt was made to force [Sigma] [less than] co when calculating the best fits, but it usually came out that way.
Initial conditions. – Luckinbill reported the exact number of organisms used to inoculate the experimental medium, so that the initial densities x(0) = 15 Paramecium/mL and y(0) = 5.833 Didinium/mL are known precisely. [[S.sub.1].sup.2], [[S.sub.2].sup.2] and [[S.sub.R].sup.2] all used these initial conditions. It was noticed, however, that much better fits could be obtained by using different initial conditions. Because uniqueness of solutions of differential equations means that the trajectories in a phase plane diagram cannot cross, in order to have diverging oscillations with successively higher peaks a trajectory must also have successively lower troughs, but the initial conditions reported for the Paramecium are as low and for the Didinium are lower than any later densities. This may reflect the fact that the Paramecium and Didinium were grown under different conditions before they were transplanted to the experimental medium, and their growth during the first half day may represent a delayed response to their previous conditions. Indeed, Luckinbill reports that he initially placed “well fed Didinium” in the experimental medium.
If the trajectories start at [t.sub.0] = 0.5 d, then the initial conditions are not known precisely, since there are undoubtedly measurement errors in the densities measured at that time, which makes x([t.sub.0]) and y([t.sub.0]) two more parameters that can be estimated. The resulting sum of squares errors, reported as [[S.sub.3].sup.2] in Table 1, are so much smaller than [[S.sub.1].sup.2] in some cases that the difference seems to be significant, not just the result of adding two more parameters. The models with time delays do not show as much difference between [[S.sub.1].sup.2] and [[S.sub.3].sup.2] because they already contain parameters that can help to fit early growth rates to the residual effects of the previous conditions. Fig. 9a, b shows that for the best model found, the least squares best fit trajectory with estimated initial conditions appears very similar to that with fixed initial conditions.
Another way to measure the error would be to compute x([t.sub.i]) and y([t.sub.i]) in Eq. 7 by resetting x([t.sub.i – 1]) and y([t.sub.i – 1]) equal to the observed values at time [t.sub.i – 1] and using these as initial conditions to solve the differential equations over only the time interval [t.sub.i – 1] to [t.sub.i], which would certainly be the correct method if there were no observation error. This “one step at a time error” was tried, but did not produce significant differences in the sum of squares errors from one model to another. Nor would the estimated parameter values produce trajectories that were anything like the observed densities if the differential equations were solved over the entire time interval without resetting x and y at each observation time.
Accuracy of parameter estimates. – The Marquardt-Levenberg method, as with any nonlinear numerical optimization procedure, searches through parameter space for the combination of parameters that gives the minimum for [S.sup.2], but cannot distinguish a local from a global minimum. The only test is to start with a different initial guess for the parameters and see if the procedure returns to the same minimum. Although it was observed that the program sometimes hung up on a local minimum, these were always obviously poor fits. A case was never observed where there were two very different sets of parameters that both gave good fits to the observed data. Near the minimum, however, small changes in some parameters may compensate for small changes in other parameters with very little difference in [S.sup.2]; the problem might be described as finding the lowest point in a narrow curving canyon that has steep sides but is fairly flat along its bottom. Therefore it seemed neither practical nor important to estimate the parameters to any high degree of accuracy, and any of them might have an error of up to 5%. The goal was not to determine the values of the parameters but what types of behavior take place in the predator-prey interaction.
Quantitative Fit of the Experimental Data: Results
A mathematical model with unspecified parameters, such as Eqs. 1-3, can be considered a hypothesis that does not predict a specific trajectory for the populations densities, but rather a family of possible trajectories, one for each set of reasonable parameter values. The hypothesis is tested by checking whether any of these trajectories closely matches the observed data. If not, the hypothesis must be rejected or modified. This process should identify models that include important features of a predator-prey interaction.
Standard model. – The best fit to the experimental data that can be obtained with the standard model, Eqs. 1-3, which incorporates a carrying capacity for the prey and a Holling Type 2 functional response, is shown in Fig. 5. Although this model correctly predicted the qualitative outcome of the experiment and produces the correct number of oscillations with reasonable sized parameters, the fit obtained is disappointing with [[S.sub.1].sup.2] = 236 137. The peaks are much too narrow and the dips are much too low and too long; it appears that both species would have died out by t = 9.
Predator mutual interference. – In Fig. 5 the predators appear to grow too rapidly when the prey is abundant. If one plots the three observed points [[x.sub.observed]([t.sub.i]), [y.sub.observed]([t.sub.i])] where the predator y(t) is at a local maximum, these points should lie on the predator isocline. Although three points are not enough to do any statistical curve fitting, it appears that this isocline is not vertical as predicted by Eq. 4b, but curves slightly to the right. Both these observations would suggest some type of predator mutual interference or intra-specific competition. This could be modeled by adding a term -[Delta][y.sup.2] to Eq. 2, to represent competition among the Didinium for space or some resource other than the Paramecium, or by replacing Eqs. 1 and 2 with
dx/dt = [Rho](1 – x/K)x – [Omega]f(x) g(y) (10) dy/dt = [Sigma]f(x) g(y) – [Gamma]y, (11)
where f(x) is still given by Eq. 3 and g(y) is given by
g(y) = y/1 – [Beta]y (12)
with [Beta] [greater than] 0, to represent a reduction in the predation rate at high predator densities due to mutual interference among the predators while searching for food. Salt (1974) also found that the predation rate per individual Didinium decreased with increasing Didinium density. Either type of predator intraspecific competition will improve the fit dramatically, but feeding interference gives a slightly better fit and seems to describe the experiment better, since there is no indication that any other resource is limiting for the Didinium.
Some authors (e.g., Freedman 1979) used g(y) = [y.sup.[Beta]] with 0 [less than] [Beta] [less than] 1 to model feeding interference, but Eq. 12 has the advantage that y/(1 + [Beta]y) [approximately equal to] y for small y, whereas [y.sup.[Beta]]/y [approaches] [infinity] as y [approaches] 0. Other authors (DeAngelis et al. 1975, Hassel et al. 1976) have used a functional response that depends on both x and y,
f(x, y) = x/[Phi] + [Beta]y + x.
This form was also tried and did not give significantly different results from Eqs. 10-12.
The best fit obtained using Eqs. 10-12 and 3 is shown in Fig. 6 and gives [[S.sub.1].sup.2] = 120 313. Although this is a great improvement over the standard model, there is still considerable systematic error. The predator peaks occur too soon and the population levels drop too low.
As a further test, reducing [Phi] to [Phi] = 3.0 gives a fair fit to the data of Fig. 1, where the predator searching efficiency was high because there was no methyl cellulose. Also keeping all parameters the same but reducing K to K = 400, to match the reduced carrying capacity shown in Fig. 3 when the Paramecium are grown alone in the half-strength cerophyl, gives a fair fit to the densities of Paramecium and Didinium grown together in Fig. 3 [ILLUSTRATION FOR FIGURES 6b, c OMITTED]. Thus on a quantitative level this model matches the results of the modifications in the experimental medium, just as the theory says it should.
Ratio-dependent functional response. – Currently there is considerable interest in the notion that the predator functional response does not depend on prey density x alone, but on the ratio of prey to predator x/y. (See the Special Feature on ratio-dependent predator-prey theory, Berryman and Matson, editors, 1992, especially the articles by Arditi and Saiah 1992, Berryman 1992, and Ginsburg and Akcakaya 1992.) Using a functional response of the form
[Mathematical Expression Omitted]
in Eqs. 1-2 gives a model that is much too stable to represent Luckinbill’s system. Although it is possible to find parameter values that produce oscillations, they do not fit the data well and the best fit to the data of Fig. 2 only had two peaks instead of three. A ratio-dependent functional response definitely does not occur in Luckinbill’s experimental system. This does not disagree, however, with the hypothesis of Arditi and Saiah (1992) that ratio dependence is caused by spatial heterogeneity.
Leslie-type model. – Leslie (1948) proposed a predator-prey model where the predator has a carrying capacity that is proportional to the prey density. Eq. 2 is replaced by
dy/dt = [Sigma](1 – [Gamma]y/x)y. (14)
This model (Eqs. 1.3. and 14) gave a better fit than the standard model with [[S.sub.1].sup.2] = 143 491, but when feeding interference among the predators was added (Eqs. 3, 10, 12 and 14), the fit was only improved slightly to [[S.sub.1].sup.2] = 129 788, with two of the error measurements in Table 1 being slightly above those of the standard model with feeding interference and two slightly below. The Leslie model with predator interference was very sensitive to the parameter [Gamma]; a 7% reduction in [Gamma] produced damped oscillations. Neither model eliminated the systematic errors noted above that the peaks are too narrow and the troughs are too low. Nor does a Leslie-type model predict the results of adding methyl cellulose or reducing the nutrients as does the standard model with predator interference. K = 400 in the Leslie-type model with predator interference gives a stable equilibrium at very low population densities.
Sigmoid functional response. – In all the variations above the models predict that the prey densities drop to much lower values than seen in the observed data. This would suggest that the Didinium reduce their searching rate when the Paramecium become scarce enough, i.e.. that the Didinium have a sigmoid (or Holling Type 3) functional response. The usual way to represent a sigmoid response is to use
f(x) = [x.sup.[Alpha]]/[[Phi].sup.[Alpha]] + [x.sup.[Alpha]] (15)
with [Alpha] [greater than] 1. (Most popular is [Alpha] = 2.) This function produces a very symmetric graph with the bottom half of the curve a mirror image of the top half of the curve. A more flexible sigmoid function, which can curve up quickly at the bottom but approach the top horizontal asymptote more slowly, is given by
f(x) = x/[Phi] + x[1 – (1 + vx)[e.sup.-vx]]. (16)
The model with Eqs. 10-12 with f(x) defined by Eq. 16 fits the data with [[S.sub.1].sup.2] = 93 546, but interestingly the best fit is with [Beta] = 0.00057 and is only slightly better than the best fit with [Beta] = 0, i.e., without any predator mutual interference. It is not clear just why a sigmoid functional response makes predator mutual interference unnecessary or even undesirable in the model, but whenever these two features were used together the fitting procedure reduced [Beta] to 0 or to a very low level. Probably, since both of these features are stabilizing, both of them together damp the oscillations too much.
Although [S.sup.2] was smaller with a sigmoid functional response than with predator mutual interference, visually the fit is not as good, at least for the predator densities [ILLUSTRATION FOR FIGURE 7 OMITTED]. Also, reducing [Phi] will not produce a good fit to the observed density curves when there is no methyl cellulose [ILLUSTRATION FOR FIGURE 7b OMITTED]. The sigmoid functional response makes the predation rate too low at low prey densities to drive the Paramecium to extinction, but even setting v = 2.0, which makes the functional response hardly different from Eq. 3, the model will still not fit the data of Fig. 1.
Delayed predator numerical response. – Comparison of the output of any of these models with the experimental data shows that the predator density in the models responds more quickly to changes in the prey density than the real Didinium do, making the peaks in graphs of predator density occur too early and be too narrow. This suggests a delay or time lag in the predator numerical response. Time lags, which are often suggested as the cause for unexplained behavior but seldom analyzed, are known to generally have a destabilizing influence and could help cause the diverging oscillations. Rather than just insert an ad hoc time lag, a mechanism will be proposed and modeled that produces the delay in the response.
Instead of assuming that the predator reproduction rate responds instantly to prey consumption, assume that (1) prey consumption creates an inflow of energy and/or nutrients into the total energy and/or nutrients stored in the predators, and (2) the overall predator-specific growth rate is an increasing function of the average energy or nutrients per predator. The simplest way to model these assumptions is to let Z be the total energy/nutrient storage of the predators and form a system of three differential equations:
dx/dt = [Rho](1 – x/K)x – [Omega]f(x) g(y), (17)
dZ/dt = [Sigma]f(x) g(y) – [Delta]Z, (18)
dy/dt = ([Lambda]Z/y – [Gamma])y. (19)
The term [Delta]Z in Eq. 18 represents the rate of energy expenditure, Z/y in Eq. 19 is the average energy/nutrients per predator, [Gamma] is, as before, the specific death rate of the predators, and [Sigma] and [Lambda] are proportionality constants with labels energy per prey and predators per energy per time, respectively. The functions f and g are defined as above. Rescaling the equations by substituting z = [Lambda]Z simplifies the equations to
dx/dt = [Rho](1 – x/K)x – [Omega]f(x) g(y), (20)
[Mathematical Expression Omitted],
dy/dt = z – [Gamma]y, (22)
where [Mathematical Expression Omitted]. Eqs. 20-22 have the advantage that there is one less parameter to estimate than in the previous set, but still require two new parameters to be estimated, [Delta] and z(0).
The variable z has the label predators per time and represents the predator’s reproductive rate. Eq. 21 says that the predator’s reproductive rate is not instantly a function of the rate of energy/nutrient intake but the rate of change of the reproduction rate is a function of the rate of energy/nutrient intake. For a given x and y the equilibrium level of z is given by
[Mathematical Expression Omitted]
and since Eq. 21 could be written dz/dt = [Delta]([z.sup.*] – z), the larger value of [Delta], the faster z will move toward [z.sup.*]. Of course [z.sup.*] is not constant since it depends on x and y. Note also that if we replaced z by [z.sup.*] in Eq. 22 and identified [Sigma] with [Mathematical Expression Omitted], Eqs. 20 and 22 would be equivalent to Eqs. 10 and 11. Thus, when [Delta] is large, z will track [z.sup.*] closely and the resulting trajectory will be very close to that given by Eqs. 10 and 11 with very little time delay in the numerical response, but if [Delta] is small, z will lag behind [z.sup.*], producing a significant delay in the predator’s numerical response to the prey density.
The model consisting of Eqs. 20 through 22 with f(x) defined by Eq. 3 and g(y) = y fits the data with [[S.sub.1].sup.2] = 72 436, but the model with g(y) defined by Eq. 12 (predator mutual interference) fits the data with [[S.sub.1].sup.2] = 51 717 [ILLUSTRATION FOR FIGURE 8 OMITTED]. However, the best fit found was Eqs. 20 through 22 with f(x) given by Eq. 16 (sigmoid functional response) and g(y) = y, which fits the data with [[S.sub.1].sup.2] = 30 084 [ILLUSTRATION FOR FIGURE 9 OMITTED]. Fig. 9b, c shows that the best fit trajectory for this model appears the same no matter which error criterion is used. Both of the latter two models give a reasonable fit to the observed data for reduced nutrient supply if K is reduced to K = 400 [ILLUSTRATION FOR FIGURES 8c AND 9e OMITTED], but to make the model with delayed numerical response and sigmoid functional response come at all close to fitting the data without methyl cellulose [ILLUSTRATION FOR FIGURE 9d OMITTED], it was necessary to not only set [Phi] = 3.0 but also [Nu] = 0.5, which makes the functional response only slightly sigmoid. This indicates that thickening the medium with methyl cellulose not only decreases the searching efficiency but also makes the Didinium reduce their search rate at low prey densities even more when the effort of swimming is greater, thus making the functional response more sigmoid than in the unthickened medium.
Comparing the values of [[S.sub.3].sup.2] (least squares error when [t.sub.0] = 0.5 and the initial conditions are estimated) in Table 1, one does not see as dramatic an improvement from the addition of the delayed numerical response. It appears that the Didinium transplanted to the experimental medium had a high average energy storage and reproduced rapidly during the first half day, which the models do not have to match when using the error criterion [[S.sub.3].sup.2]. Still, the delayed numerical response does more than improve the fit at the early times, since adding it to the sigmoid functional response reduces [[S.sub.3].sup.2] from 41 775 to 26 616 and Fig. 9 is clearly a better fit overall than Fig. 7.
The model with delayed numerical response and sigmoid functional response was the best one found based on the criteria of best fit with the least complexity. Adding predator mutual interference to the model with delayed numerical response and sigmoid functional response gave [[S.sub.1].sup.2] = 29 231, which is not enough improvement to justify the additional parameter to he estimated. The data for the Paramecium grown alone [ILLUSTRATION FOR FIGURES 2 AND 3 OMITTED] show that the densities oscillate around the carrying capacities, suggesting that there might also be delays in the prey growth rate, but a model that also introduced a delay in the prey growth rate in a similar manner and used four equations with 10 parameters did not improve the fit enough ([[S.sub.1].sup.2] = 25 439) to justify the added complexity. There may be a delay in the response of the Paramecium growth rate to the food supply, but when grown together with the Didinium their density never approached a level where the food supply was limiting and hence the delay does not show up.
Discussion and Conclusions
A standard predator-prey model with a carrying capacity for the prey and a saturating functional response gives a qualitative explanation of the results of Luckinbill’s (1973) experiment that both thickening the medium and reducing the nutrient supply help stabilize the oscillations, but gives only a poor quantitative fit to the data. The quantitative fit can be dramatically improved by modifying the model to include predator mutual interference or a sigmoid functional response, but strangely not both. The best fit, however, is achieved by coupling either of these with a predator growth rate that is proportional to nutrient and/or energy storage, which in turn has a rate of change determined by prey consumption, thus producing a delayed numerical response (see Table 1 and [ILLUSTRATION FOR FIGURES 8 AND 9 OMITTED]). Other modifications either gave poor fits, were overly sensitive to the parameter values, or required unreasonable parameter values.
When populations are measured in terms of biomass, then it makes sense that the population growth rate would depend directly on the rate of nutrient intake, but when populations are measured in terms of numbers of individuals as in Luckinbill’s (1973) study, then it makes sense that the reproductive rate would depend on the nutrient storage and thus in a delayed manner on the nutrient intake. This characteristic was certainly important for modeling the Didinium and may be important for many other species. It deserves further study.
The models produced by coupling the delayed numerical response with predator mutual interference or with a sigmoid functional response both give excellent fits, but the latter seems to be the best. The former, after estimating parameters to fit the data of Fig. 2, also gives a good fit to both the data without methyl cellulose and with reduced nutrients by making the expected changes in [Phi] and K [ILLUSTRATION FOR FIGURE 8 OMITTED]. Some of the parameter values for this model seem unusually large since they are higher than found in any of the other models, but [Rho] = 3.40 is not too much higher than the a priori estimate and [Omega] = 43.37 matches the maximum predation rate reported by Salt (1974). In contrast, the parameter values for the latter model all seem very reasonable, with the value [Omega] = 9.74 being close to Reukauf’s (1930) value of 12 for the maximum predation rate, but to accept it as the best model also implies accepting that the functional response becomes more sigmoid as the medium is thickened [ILLUSTRATION FOR FIGURE 9 OMITTED]. It is not unreasonable to expect that in the thicker medium the Didinium quit or at least reduce their searching for prey when the prey densities are low, but that has not been verified independently.
Overall, the ability to fit actual experimental data for a two-species interaction with an [R.sup.2] of 0.94 using only a few differential equations is remarkable and validates the utility of these models. The fact that after fitting the parameters to one set of data the model could predict quite well the results of two other experiments further strengthens our confidence in the model. The results seem to depend more on the features being modeled than on the specific mathematical formulations used. For example, four different mathematical formulas to represent predator mutual interference were tried and all gave about the same results.
Of course, some systematic error still remains. All of the models still predict that the population levels of the Paramecium drop too low between the peaks, and none of them clearly predicts that the Didinium will go extinct before the Paramecium, as was always observed by Luckinbill. It would be instructive to make a size-structured model, especially since Luckinbill reports that much of the reduction in the growth rate of the Didinium when the nutrients for the Paramecium were reduced is attributable to the smaller average size of the Paramecium, but there does not seem to be enough information on how reproduction varies with size to do so at present.
There are also technical difficulties with using the method of least squares to fit oscillating data. The trajectories are sensitive to initial conditions, as shown by the remarkably better fits attainable by starting at [t.sub.0] = 0.5 and estimating initial conditions instead of starting at [t.sub.0] = 0 and using the initial conditions reported by Luckinbill (1973) (Table 1). Small changes in initial conditions or parameters often produce trajectories that have oscillations with about the same period and magnitude, but are slightly out of phase. The graphs may appear quite similar to us, but are measured as very different by the least squares method because it only measures vertical distances in the graphs. The long data set [ILLUSTRATION FOR FIGURE 3 OMITTED] is especially hard to fit using least squares techniques because there are several places where the Paramecium density makes a brief upward jog, probably caused by the fact that Luckinbill (1973) transferred the populations to fresh medium every 1.5 or 2 d. These brief perturbations don’t change the overall pattern but do lengthen the cycles where they occur, and there is no way that the even cycles from the theoretical model can stay in phase with the nonuniform cycles of the data.
Predator-prey models such as those used in this paper are often criticized for being too simplistic and for using phenomenological parameters that cannot be measured directly. The goal of modeling is to capture the essence of the interaction without having to represent all the complexity. Two or three simple differential equations that can capture the essence of the Paramecium-Didinium interaction and explain the results of different experimental manipulations are a better model than a large complex simulation that can do the same thing. It is true that the parameters are phenomenological in the sense that they represent population characteristics rather than individual characteristics and can only be measured indirectly. We cannot measure the intrinsic growth or the predation rate of an individual Didinium; they must be calculated from changes in the populations. But the same can be said for many useful concepts in science. We cannot measure velocity directly but compute it from two measurements of position and one measurement of the elapsed time, and computing acceleration is even more complicated. Is it any worse that we have to measure predation rate from two measurements of prey population and one measurement of elapsed time? Energy is an extremely useful concept in many sciences, including ecology, but may only be a construct of our own minds to explain certain regularities we observe in natural phenomena. Nor can even such a universally accepted phenomenon as temperature be measured directly, but instead is inferred from changes in the volume of a liquid or metal. Instead of discarding these quantities as phenomenological we give them names and work at finding better ways of measuring them indirectly, because they are useful. The same should be true of the quantities represented by the parameters in a predator-prey model.
Peters (1991) has criticized much of theoretical ecology, including Lotka-Volterra models, as being nonpredictive and hence untestable. All too often this criticism is valid, and efforts at model validation are far too rare. But the type of models studied in this paper do make predictions and are testable. Even without the parameter values specified, Eqs. 1-3 predict a collection of possible trajectories. If the qualitative outcome of an experiment does not match any of these possible trajectories, then the model must be discarded or overhauled. Thus the original predator-prey model of Lotka and Volterra, which predicts only neutrally stable cycles under all conditions, has long since been discarded and replaced by an improved version, Eqs. 1-3, which predicts shifts between rapid extinction, diverging oscillations, and asymptotically stable limit cycles as the experimental setup is modified. Another prediction, that further thickening of the medium or reduction of the nutrients would lead to stable equilibrium populations, has yet to be tested. The notion of a least squares best fit is to find the trajectory out of the predicted collection of possible trajectories that comes closest to fitting the observed data. If this best fit trajectory does not quantitatively match the data, the others will match even worse and the model must again be rejected or modified. On this level Eqs. 1-3 did not fare too well, and various modifications were tried. A model with sigmoid functional response and delayed numerical response due to the predator reproduction rate varying with energy/nutrient storage and not with rate of prey consumption is the best model to date. Since this is not only a general model, but also has parameter values estimated, it makes very specific predictions. It correctly predicted the outcome of the low nutrient experiment. An unexpected prediction from comparing this model to the unthickened medium experiment is that the parameters in a sigmoid functional response depend on the relative cost of searching to the benefit of capturing prey. As searching becomes more difficult the functional response becomes more sigmoid, that is, the search rate at low prey densities is more depressed. Experiments to test this hypothesis directly could greatly increase our understanding of the predator functional response. One of the key values of modeling is the identification of such opportunities for further experiments.
Further testing should be done, but the models developed seem adequate for the Paramecium-Didinium interaction in a simple environment. As this and similar models are tested for other species and other environments, they will probably need additional modifications, but they have promise to become a valuable theory for explaining and predicting population levels in predator-prey interactions.
Acknowledgments
The author wishes to thank Don DeAngelis, Carl Walters, and an anonymous reviewer for helpful comments on the manuscript, and the Department of Mathematics at the University of Tennessee for partial support while working on this project.
Literature Cited
Arditi, R., and H. Saiah. 1992. Empirical evidence of the role of heterogeneity in ratio-dependent consumption. Ecology 73:1544-1551.
Beddington, J. R., M. P. Hassell, and J. H. Lawton. 1976. The components of arthropod predation. II. The predator rate of increase. Journal of Animal Ecology 45:165-186.
Berryman, A. A. 1992. The origins and evolution of predator-prey theory. Ecology 73:1530-1535.
Berryman, A. A., and P. Matson, editors. 1992. Special feature: ratio-dependent predator-prey theory, Ecology 73: 1529-1566.
Coe, M. J., D. H. Cumming, and J. Phillipson. 1976. Biomass and production of large African herbivores in relation to rainfall and primary production. Oecologia (Berlin) 22: 341-354.
Danby, J. M. A. 1985. Computing applications to differential equations. Reston Publishing, Reston, Virginia, USA.
DeAngelis, D. L., R. A. Goldstein, and R. V. O’Neill. 1975. A model for trophic interaction. Ecology 56:881-892.
Flanders, S. E. 1948. A host-parasite community to demonstrate balance. Ecology 29:123.
Freedman, H. I. 1976. Graphical stability, enrichment, and pest control by a natural enemy. Mathematical Biosciences 31:207-225.
—–. 1979. Stability analysis of a predator-prey system with mutual interference and density-dependent death rates. Bulletin of Mathematical Biology 41:67-78.
Gause, G. F. 1934. The struggle for existence. Williams and Wilkins, Baltimore, Maryland, USA.
Gause, G. F., N. P. Smaragdova, and A. A. Witt. 1936. Further studies of interaction between predator and prey. Journal of Animal Ecology 5:1-18.
Ginzburg, L. R., and H. R. Akcakaya. 1992. Consequences of ratio-dependent predation for steady-state properties of ecosystems. Ecology 73:1536-1543.
Harrison, G. W. 1986. Multiple stable equilibria in a predator-prey system. Bulletin of Mathematical Biology 48: 137-148.
Hassell, M. P., J. H. Lawton, and J. R. Beddington. 1976. The components of arthropod predation. I. The prey death rate. Journal of Animal Ecology 45:135-161.
Hastings, A. 1978. Global stability of a two species system. Journal of Mathematical Biology 5:399-403.
Holling, C. S. 1959. The components of predation as revealed by a study of small-mammal predation of the European pine sawfly. Canadian Entomologist 91:293-320.
Hsu, S. B. 1978. The application of the Poincare transform to the Lotka-Volterra model. Journal of Mathematical Biology 6:67-73.
Huffaker, C. B. 1958. Experimental studies on predation: dispersion factors and predator-prey oscillations. Hilgardia 27:343-383.
Leslie, P. H. 1948. Some further notes on the use of matrices in population mathematics. Biometrica 35:213-245.
Lotka, A. J. 1925. Elements of physical biology. Williams and Wilkins, Baltimore, Maryland, USA.
Luckinbill, L. S. 1973. Coexistence in laboratory populations of Paramecium aurelia and its predator Didinium nasutum. Ecology 54:1320-1327.
Murdoch, W. W., and A. Oaten. 1975. Predation and population stability. Advances in Ecological Research 9:1-131.
Peters, R. H. 1991. A critique for ecology, Cambridge University Press, Cambridge, England.
Press, W. H., B. P. Flannery, S. A. Teukolsky, and W. T. Vetterling. 1986. Numerical recipes. Cambridge University Press, Cambridge, England.
Reukauf, E. 1930. Zur biologie von Didinium nasutum. Zeitschrift fur vergleichende Physiologie 11:689-701.
Rosenzwieg, M. L. 1969. Paradox of enrichment: destabilization of exploitation systems in ecological time. Science 171:385-387.
Rosenzwieg, M. L., and R. H. MacArthur. 1963. Graphical representation and stability conditions of predator-prey interactions. American Naturalist 47:209-223.
Salt, G. W. 1974. Predator and prey densities as controls on the rate of capture by the predator Didinium nasutum. Ecology 55:434-439.
Slobodkin, L. B. 1986. The role of minimalism in art and science. American Naturalist 127:257-265.
Volterra, V. 1931. Variations and fluctuations of the number of individuals in animal species living together. Translated from 1928 edition by R. N. Chapman. Animal ecology. Arno, New York, New York, USA.
Wichterman, R. 1986. The biology of paramecium. Plenum, New York, New York, USA.
COPYRIGHT 1995 Ecological Society of America
COPYRIGHT 2004 Gale Group
