May 20, 2018, Sunday
University of Colorado at Boulder Search A to Z Campus Map University of Colorado at Boulder CU 
Search Links

Extensions to The Dynamics of Arthropod Predator-Prey Systems

From MathBio

Jump to: navigation, search

Review by Justin Baacke, Nicole Look and Brit Schneiders, May 2012.

The following wiki reproduces the results presented by one of the insect host-parasite models developed in Arthropod Predator-Prey Systems by M. P. Hassell. More specifically, the article below first demonstrates the analytical and computational analysis of the model for non-random search and constant searching efficiency. The article then extends the model by in two ways: by extending from a discretized model to continuous space of non-random search, and also incorporating hyperparasitoidism. The work shown below can be attributed to Hassell unless otherwise specified.

Executive Summary

The non-random search model presented below, as well as other models that describe species interactions, developed from modifying the Lotka-Volterra predator-prey model of 1910. Although the Lotka-Volterra model simplistically describes the interactions between two species, it only considers continuous time. Nicholson and Bailey (1935) resolved this problem by creating a predator-prey model for discrete time. More recently, Hassell and May proposed the replacement of predators with the notion of parasitoids to create simpler models with new extensions. A parasitoid can be thought of as the cross between a parasite and a predator: it is like a parasite in physiology, but like a predator in the ecological outcome of killing its host, as a predator kills its prey. While parasitoids span a variety of phylums, Hassell specifically investigates the Arthropoda phylum.

Beyond using parasitoids to simplify their model, Hassell and May also looked at a few different extensions, one of which was Non-Random search. Random search is an easier model to create, but non-random search is a logical modification because it is seen all around in nature. In the following article, we reproduce both the random and non-random search models in discrete time. Then we also use an agent-based model to simulate what might happen for both random and non-random search in continuous time. Lastly, we introduce the concept of hyperparasitoidism, and look at the effects of this concept on our models. Hyperparasitoids are parasitoids that infect other parasitoids, so essentially we are changing our two-species system in a three-species system in which there are two levels of predators, and two levels of prey. We found both our random search models, discretized and continuous time, to be unstable for both two-species and three-species systems. Then, when we implemented non-random search, we found stable results in discretized and continuous time, for both parasitoids vs. hosts, and hyperparasitoids vs. parasitoids vs. hosts.

Biological Context

Mathematically, parasitoids are a good, simple model of predators when describing predator-prey interactions, because, like predators, they kill their host and also because they remain in confined spaces. Parasitoids of arthropods (arachnids, insects, crustaceans, etc) are especially great models because in many species, one parasitoid produces one offspring per every host it kills. This one-to-one ratio between the average number of predator progeny produced and prey attacked becomes important in the upcoming model. Simplifications such as this make the upcoming model easy to apply.

Even within Arthropoda, there is an extensive diversity of parasitoids, as can be seen in the figure below:

Figure 1: Arthropod Parasitoids Taxonomy

Figure 1: Within Arthropoda, there are many orders that include parasitoids, the most common of which is Hymenoptera: namely, wasps and bees. An estimated 74% of Hymenapterans are parasitoids, and as an order, hymenapterans make up about 75% of all discovered parasitoids. Next most common is the order Diptera, where you find common flies. The third order of Arthropoda that includes many parasitoids is Coleoptera, comprised of many species of beetles.

As mentioned, a parasitoid represents a cross between a parasite and a predator. By definition, a parasite is an organism that lives in or on another organism (host) and benefits by obtaining nutrients from the host at the host’s expense. Because the parasite relies on the host for survival, the parasite generally avoids killing the host. In contrast, a predator survives by preying on other organisms (prey) that result in their death. Therefore, a parasitoid mimics a parasite during its larval stage, when it feeds on the body of the host, but represents a predator since the parasitoid eventually kills the host.

Most parasitoids lead protelean life cycles, meaning they are parasitic as larvae and are free living as adults. While each parasitoid’s life cycle may vary, it generally follows as: the free-living female adult releases its eggs into the environment to be ingested or injects the eggs directly into the next host; the eggs develop into larvae while feeding on the host; the larvae develop into juveniles, exit and kill the host. The picture below demonstrates a typical hymenopteran life cycle (ants, bees, wasps):

Figure 2: Hymenapteran Life Cycle

Figure 2: Starting at the top left, an adult female oviposits her egg inside of the host (generally an arthropod), her egg hatches and gives rise to a larva, using much of its host's nutrients, the larva develops into a pupa, usually taking up most of the host's body cavity. The pupa develops into a young adult. When the adult is ready, it emerges from its host, killing the host in the process. Note that this example has a one to one ratio of predator progeny produced per prey attacked.

Parasitoids that exhibit this one to one ratio are called "Solitary Parasitoids". Most hymenopteran parasitoids are solitary as opposed to "Gregarian", meaning the parasitoid produces more than one progeny per host parasitzed, or "Polyembryonic", meaning one egg is oviposited, but the egg splits into many larva. Below are examples of each:

Figure 3: A gregarious parasitoid: generally an ectoparasitoid (meaning the larvae develop externally of the host) that produces more than one larva per host infected.
Figure 4: A polyembryonic parasitoid: one egg becomes many larvae.
Figure 5: A solitary parasitoid: one larva develops inside the host.

Another biological concept that we will explore is that of hyperparasitoids. Hyperparasitoids, in terms of a predator-prey system food web, are analogous to "tertiary consumers". A primary consumer is generally an herbivore, and the first prey. A secondary consumer predates that prey, and in turn serves as prey to the tertiary consumer. Think of a hawk, which eats a snake, which preys on a mouse (which consumes plants like grass).

In the same sense, a hyperparasitoid (or secondary parasitoid) parasitizes a primary parasitoid, which in turn parasitizes its host. This is an important element of parasitoid systems, and was not well-studied at the time Hassell publicized his book on Arthropod Predator-Prey Systems.


Predator-prey interactions were originally explained by the Lotka-Volterra model, proposed in 1910. However, many researchers found it necessary to redefine predator behavior assumptions in order to account for various forms of interaction. In 1935, Nicholson and Bailey modified the Lotka-Volterra model to incorporate discrete time. More specifically, the model describes the population dynamics of a coupled host-parasite system. While it assumes that parasites search hosts at random and that the environment contains an even distribution of hosts and parasites, it does not allow for stable host-parasite interactions. In addition, it yields diverging oscillations after perturbations. Most species populations do not tend toward extinction. For a more in-depth explanation of the Nicholson-Bailey model, please consult “The Dynamics of Arthropod Predator-Prey Systems.”

To add stability and introduce non-random search, Hassell and May (1973) then modified the Nicholson-Bailey model. By considering a nonhomogeneous distribution of prey and a non-random search path by the predators, Hassell and May reconciled the stability problem. The modified model’s assumptions are reasonable for two main reasons. First, numerous examples of such behavior have been observed in nature. For example, ladybugs (predator) tend to remain on leaves with aggregated populations of aphids (prey) rather than leaves without aphids. Second, laboratory experiments using these assumptions more frequently reproduce stable populations than experiments that include homogenous spatial prey distribution.

Mathematical Model: Reproducing Results

The majority of insect host-parasite models follow a general form:

                                           P_{{(t+1)}}=N_{{t}}-N_{{s}}                    (1)

where N_{{s}} is the number of host survivors once P_{{t}} parasites have searched for N_{{t}} hosts that results in a subsequent generation of P_{{(t+1)}} parasite progeny. Modifications of the general model can be made by changing the function f[P_{{t}},N_{{t}}], which also contains all assumptions made about the parasite searching behavior. Hassell and May defend models with fewer parameters and fewer assumptions as more applicable to natural populations. They further propose that the stability of most natural host-parasite interactions is affected by non-random search. This proposition contradicted the previous host-parasite models, including that of Nicholson and Bailey (1935), and led to an aggregative model.

The aggregative model assumes that certain areas contain higher densities of prey than other areas. Therefore, it becomes more profitable for predators to develop a non-random search pattern and take advantage of the high-density patches. Fruit flies exemplify this behavior, such that their feeding sites are more likely to be clumped together. Investigators found that hungry flies would change their behavior by taking larger steps and smaller turns to cover a greater area. Shortly after finding food, step length then shortened and average angle turn increased so that the fly could explore the area with greater detail where the food had been found. Refer to the Example section of The Dynamics of Arthropod Predator-Prey Systems wiki for more explanation.

Hassell and May’s model incorporates ideas of aggregation and therefore provides far greater stability. The model considers "n" different area patches, which may or may not contain a number of prey. Each patch is not necessarily a predefined geometric spatial region since the patch must be linked with the modified foraging path of the parasitoid. The patch modification is presented by the following function that can then be incorporated into the general system of equations (1)

                                      f\left(N_{t},P_{t}\right)=\sum _{{i=1}}^{n}\alpha _{i}e^{{-a\beta _{i}P_{t}}}          (2)

with the newly defined variables:

                                      \alpha _{i}    Fraction of host population in patch i
                                      \beta _{i}    Fraction of parasite population in patch i 
                                      n     Number of patches
                                      F     Prey rate of increase

Equation (1) can then be modified as:

                                       N_{{s}}=N_{{t}}\sum _{{i=1}}^{n}{[\alpha _{{i}}e^{{-a\beta _{{i}}P_{{t}}}}]}        (3)

The summation of the considered space allows for prey to non-uniformly occupy any of the patches, such that P_{{t}} parasites and N_{{t}} hosts can be distributed into n areas specified by \alpha _{{i}} and \beta _{{i}}. Equation (3) assumes that within each area (i) the random search for hosts is random and the parasites’ ability to find hosts is independent of host and parasite density. Furthermore, it is important to note that when the predators are evenly distributed, the \beta _{i}={\frac  {\beta _{t}}{n}} and function (2) simplifies to the Nicholson-Bailey model.

Analytical Analysis

Since the purpose of each model is to predict the expected outcome for host and parasite populations after interaction, it is necessary to study the stability in regards to each parameter. This section will first study the stability of the general system of difference equations and then go on to incorporate function (2). When the parasite population is specified and synchronized temporally with the host population, the generalized model (1) can be written as

                                         P_{{(t+1)}}=N_{t}-{\frac  {N_{{(t+1)}}}{F}}          (4)

where N and P are the host and parasite densities in generations t and t + 1 and F is the rate of increase of the host population.

The equilibrium solution can be found by letting N_{{(t+1)}}=N_{{t}}=N^{{*}} and P_{{(t+1)}}=P_{{t}}=P^{{*}}:

                                         N^{{*}}=FN^{{*}}f[P^{{*}},N^{{*}}]\Rightarrow f[P^{{*}},N^{{*}}]={\frac  {1}{F}}
                                         P^{{*}}=N^{{*}}-{\frac  {N^{{*}}}{F}}\Rightarrow FP^{{*}}=(F-1)N^{{*}}    (5)

provided that F>1. The equilibrium solution is only meaningful if the system actually returns to its equilibrium populations after perturbed. In order to assess the stability of the equilibrium solution (N^{{*}},P^{{*}}), the perturbed populations can be written as

                                         P_{{t}}=P^{{*}}(1+y_{t})    (6)

where x_{{t}} and y_{{t}} are small initial perturbations. Given the total number of prey N_{{a}} attacked by P_{{t}} predators, the dynamics of the perturbations can be studied by Taylor expansion of

                                         N_{{a}}=N_{{t}}[1-f(N_{{t}},P_{{t}})]     (7)

about the equilibrium. After discarding higher order terms

                                         x_{{(t+1)}}=(1-\nu )x_{{t}}-\eta (F-1)y_{{t}}
                                         (F-1)y_{{(t+1)}}=Fx_{{t}}-x_{{(t+1)}}      (8)

where \eta \equiv -N^{{*}}({\frac  {\partial f}{\partial P}})^{{*}} and \nu \equiv FN^{{*}}({\frac  {\partial f}{\partial N}})^{{*}} and the partial derivatives are evaluated at the equilibrium point. Since the fraction of un-parasitized hosts most likely decreases when the number of parasites increases and increases when the number of hosts increases, it is expected that \eta >0,\nu >0.

Standard methods for analyzing stability transform the solution of the linear difference equations into the form of

                                         x_{{t}}=A_{{1}}(\lambda _{{1}})^{{t}}+A_{{2}}(\lambda _{{2}})^{{t}}
                                         y_{{t}}=B_{{1}}(\lambda _{{1}})^{{t}}+B_{{2}}(\lambda _{{2}})^{{t}}      (9)

By substitution of equations (9) into equations (8), the expanded equations can be rewritten as

                                         (1+\nu -\lambda )A-\eta (F-1)B=0
                                         (F-\lambda )A-(F-1)\lambda B=0               (10)

Equations (10) are only consistent when

                                         det{\begin{bmatrix}(1+\nu -\lambda )&(F-\lambda )\\-\eta (F-1)&-(F-1)\lambda )\end{bmatrix}}=0

leading to the characteristic equation

                                         \lambda ^{{2}}-\lambda (1+\nu +\eta )+F\eta =0             (11)

which can be solved as

                                        \lambda _{{1,2}}={\frac  {(1+\nu +\eta )\pm {\sqrt  {(1+\nu +\eta )^{{2}}-4F\eta }}}{2}}       (12). 

The overall stability criterion is |\lambda |<1:

  1. If the factor under the square root is positive, solutions \lambda _{{1,2}} are real numbers and the damping is exponential.
  2. If the factor is negative, \lambda _{{1,2}} are a conjugate pair of complex numbers and the oscillatory behavior has decreasing amplitudes.
  3. If |\lambda |>1, the solutions are either purely exponential growth, when the factor is positive, or growing oscillations, when the factor is negative.

Applying the stability criterion to solution (12) leads to the overall stability criterion of \eta , such that

                                       {\frac  {(1+\nu +\eta )\pm {\sqrt  {(1+\nu +\eta )^{{2}}-4F\eta }}}{2}}>1\Rightarrow 
                                       {\frac  {\nu }{(F-1)}}<\eta <{\frac  {1}{F}} and \eta >{\frac  {-(2+\nu )}{(F+1)}}       (13). 

It can be noted that

  1. Larger values of \eta outside of the range lead to unstable oscillations
  2. Smaller values of \eta outside of the range lead to unstable monotonic growth
  3. Within the range, stability can be monotonic when \eta <\eta _{{0}} and oscillatory when \eta >\eta _{{0}}, where \eta _{{0}}=({\sqrt  {F}}-{\sqrt  {F-1-\nu }})^{{2}}

More specific for the non-random search model, where f(N,P)=\sum _{{i}}\alpha _{{i}}e^{{-a\beta _{{i}}P}}. In this case, the equilibrium populations N^{{*}} and P^{{*}} found in equation (5) can be applied to equation (8) to yield \nu and \eta as

                                      \nu \equiv FN^{{*}}{\frac  {\partial f}{\partial N}}^{{*}}=0
                                      \eta \equiv -N^{{*}}{\frac  {\partial f}{\partial P}}^{{*}}={\frac  {F}{(F-1)}}\sum _{{i}}\alpha _{{i}}(a\beta _{{i}}P^{{*}})e^{{-a\beta _{{i}}P^{{*}}}}     (14).

The overall stability criterion (13) leads to the exact condition for stability

                                      F\sum _{{i=1}}^{n}\alpha _{{i}}a\beta _{{i}}P_{{i}}e^{{-a\beta _{{i}}P^{{*}}}}<{\frac  {F-1}{F}}            (15). 

It is important to note that stability is most affected by:

  1. F, the effective rate of increase of the host population
  2. \alpha _{{i}}, the distribution of the hosts
  3. \beta _{{i}}, the distribution of the parasites

Furthermore, the scenario analyzed by Hassell and May considered one patch with a high density fraction, \alpha , of the population and the remaining "n-1" patches containing {\frac  {1-\alpha }{1-n}} as a density fraction of the prey. This allows for the introduction of an aggregation index \mu which is defined as

                                       \beta _{i}=c\alpha _{i}{}^{{\mu }}        (16)

where c is a normalizing constant.

Here \mu is considered a measure of the distribution of predators. When \mu =0 a random search pattern exists, and when \mu =\infty , the case in which all predators are on the patch with the highest density is studied.

Hassell and May concluded that 4 elements contribute to the stability of a given model:

  1. \mu : Aggregation. It was found that higher values of predator aggregation contributed to greater stability in a given model. This finding can be used to stabilize models that are otherwise highly unstable.
  2. \lambda : Stability. Stability is highly dependent on prey rate of increase. As in all cases, models become unstable when prey rate of increase rises. Lotka and Volterra first and most notably demonstrated this in their studies of why increasing the prey fish populations led to an extinction of the long-term population.
  3. \alpha : Prey population. As more of the population concentrates on the high-density patches, higher values of \lambda allow for stability. However, values of \alpha >.5 do not permit stability at low values of \lambda . If \alpha is small and thus there is little difference between the high and lower density patches, lower values of \mu are possible for stability. Although, it is important to note that if \alpha is "too" low, then the model collapses back into the basic Nicholson model and instability resumes.
  4. n-1: The number of low prey density patches. With a fixed value of \alpha , less predator aggregation is necessary for stability because more low-density patches are available. Maynard Smith, Roff, and Hilborn (1974, 1974, and 1975, respectively) investigated this phenomenon, in which stability has a direct relation to the number of spatial sub-units.

Biological Context of Non-Random Search

Non-random search is an intuitive modification to make since that is what is seen in nature. In Arthropoda, and in Insecta specifically, adult female parasitoids are the ones who do the host finding, since the purpose of finding the host is to lay her eggs. The adult female does this in a number of ways. One way is by attraction to the host's calls. This can be seen with cicada hosts, specifically. Cicadas have very strong mating calls, which female parasitoids cue in to.

Another way females find their hosts is by volatile cues. These are chemical cues exhibited by either the host or the plant upon which the host feeds. One such example is aggregation and sex pheromones, which a species uses to attract other members of the species. Parasitoids cue in to the same chemical cue. Parasitoids also cue in to plant damage, because this leads them to the host they are searching for, due to the fact that their hosts are usually herbivorous insects. Lastly, parasitoids are attracted to the scent of frass, which is a powdery material herbivorous insects pass as waste from digesting plants. This again leads them to their host.

Computational Model

Computational Analysis

In order to reproduce the numerical results of Hassell May, we first used a discrete space numerical model in MATLAB in order to examine the behavior of a perturbed parasitoid system using both one patch (as fits the Nicholson Bailey approach) and then a multi-patch system. We were able to very clearly reproduce the results obtained by Hassell may as to both the instability of the Nicholson Bailey model and the stable return to a steady state of the Hassell May model. Below is a model of what occurs when we model a slightly perturbed parasitoid (red) host (blue) system using the Nicholson Bailey model.

Figure 6: A plot of a perturbed Nicholson Bailey Model using random search paths.

As you can see there are small oscillations towards the beginning eventually leading to an explosion of first the host population, then the parasitoid population, which causes the death of all of both species. In order to rectify this we incorporated Hassell May’s multi patch distribution into our model, when we do this we can see that our slightly perturbed steady state model oscillates towards its steady state as demonstrated in the following picture:

Figure 7: A plot of a perturbed Nicholson Bailey Model using non-random search paths.

Analysis of Convergence Time

As a sidenote for the numerically curious we also examined the convergence times when we varied various parameters of the system. The following series of plots are the result of this analysis, it is important to note that if the system did not converge for a given value then it was given a convergence time of zero. All other variables were held at a default constant during the variation of each individual parameter.

First we examined stability with respect to prey reproduction rate or F as referred to in the model. As you can see values short of 1.25 did not converge for this model as this is to be expected, and truly quick convergence was not achieved until a small range after 2, after which the time to regain stability asmytotically decreased. This reflects the stabilizing nature of higher values of F as shown in the analytical results.

Figure 8:Number of iterations to convergence based off differing values for the prey rate of increase.

Next we examined the speed of convergence with respect to mu, or the tendency of the predators to converge on the patch with the highest concentration of parasitoids. Here we see an odd behavior as for lower values of mu it converges relatively quickly until it reaches somewhere around .3 then we see a spike and a bowl of convergence inside the range of .3 to .6, this is likely the most stable range for the system as we can see the system converges for all values in this range. There is a region of non-convergence from .6 to .85 and then we resume the relatively quick convergence which is likely due to the system dying off.

Figure 9:Number of iterations to convergence based off differing values for the parasitoid aggregation constant.

The number of patches results in a curious analysis of how segmented a system needs to be in order for the Hassell May model to take effect. As we can see from the plot below, anything below around 50 plots does not result in convergence and anything until around 100 does not necessarily converge. This is due to the fact that in this range we are too close to resembling the single patch case. As we can see when we increase the number of patches we gain a more and more steady rate of convergence and do not have the instability associated with the Nicholson Bailey model.

Figure 10: Number of iterations to convergence based off differing values for the number of patches used in the model.

Alpha, or the tendency of the hosts to congregate in one patch can also be shown to have a nice range of stability for a given system. We see unstable behavior outside the range from about .3 to .6 however inside this range we have a fairly steady convergence rate, this would indicate that this is the preferable range for the system in order to return to steady state.

Figure 11:Number of iterations to convergence based off differing values for the aggregation constant for parasitoid.

Finally Area of discovery, or a in the model, was analyzed. Here we can see this system exhibits a very interesting behavior as there seem to be several different zones with stability associated with them. The most prominent (however leading to the largest time to convergence) is roughly .05 to .175 inside of which the system converges nicely. Both before and after this we see an almost periodic tendency to converge. This is indicative of drastic changes in the steady state of the system as, were we to look at the plots of these systems, we would see that they have a tendency to have much higher, or much lower steady states, depending of which side of the (.05, .175) range they are on.

Figure 12:Number of iterations to convergence based off differing values for the area of discovery.

Agent Based Models

In order to extend our model into continuous space we also analyzed the effect of a random search path, versus a non-random search path in an agent based modeling system. An agent based modeling system is one in which the modeler programs in a set of “base behaviours” for the various entities in their system to follow (based off research on the given system) and then generates a large number of these entities and observes the interaction. While this is not the most mathematically rigourous method of dealing with a system, it allows the users to observe how interactions vary when throughout large systems of entities with varying “base behaviours.”

In order to do this the agent based modeling platform Netlogo was used and the hosts and parasitoids were designed to follow path algorithms similar to those described by Hassell and May. For the random search both the hosts and the parasitoids used a random search path, if a parasitoid came in range of a host then the host would become a parasitoid. The parasitoids eventually died of age. As predicted by the discrete model the populations both eventually exploded and both died off.

In the nonrandom search model the hosts tended to congregate towards the center of the map and the parasitoids shortened their step size and widened their average turning angle for a short period after they had encountered a host. This model resulted in what could best be described as a tendency to move towards what resembled a steady state for a stochastic system. While not remaining at one value the populations of the hosts and parasitoids did not vary considerably over time.

Due to the limitations of what someone can put on a Wikipedia page these models are available upon request from

Extension: Hyperparasitoids

Hyperparasitoids, as briefly mentioned before, are parasitoids that prey on other parasitoids, and can be thought of as "tertiary consumers" in a typical food web.

A predator-prey system, or in this case, a parasitoid-host system, is oversimpified in the sense that it doesn't take into account many factors that also play into a species population. For example, the prevalence of the primary parasitoid (or secondary consumer in the predator-prey model) undoubtedly depends on the population of its prey, but also depends on the population of its respective predator. Studying hyperparasitoids within this system can shed light on how more complex species dynamics play into the prevalence of each respective species. Not only would the introduction of hyperparasitoids affect the primary parasitoid population, but also the original host population, since it is dependent on the primary parasitoid population as well.

Here in this picture, you can see a number of ectoparasitoids infecting a large, green catepillar. Circled in red is a hyperparasitoid infecting one of the ectoparasitoids!

Figure 13: A hyperparasitoid atop its parasitoid host.

Discrete Computational results for Hyperparasitoids

Here we will first look at a discrete model of a random search hyperparasitoid (green), parasitoid (red), host (blue) system as would be described by Nicholson Bailey. As you can see this model quickly fails as all of the species die off in very few time steps when they are perturbed from their steady states.

Figure 14: A plot of population over time for a hyperparasitoid, parasitoid, host system using the Nicholson Bailey mode with a random search pathl

Fortunately when we consider the model for a nonrandom search for hyperparasitoides as would have been implemented by Hassell May and is generated easily out of a few modifications to the code used above for the nonrandom search for a parasitoid system, the system stabilizes into what one might consider “oscillatory steady states” as shown below. This impressive dynamic nicely demonstrates the population dynamics for a nonrandom search for a hyperparasitoid model.

Figure 15: A plot of population over time for a hyperparasitoid, parasitoid, host system using the Nicholson Bailey mode with a non-random search path

Agent Based Hyperparasitoid Model

In order to continue on our extension to continuous space we also modeled hyperparasitoides in an agent based modeling scenario for both random and non-random search patterns. There are however a few differences notable from the parasitoid host model. In addition to the hyperparasitoides, two “holding tank” species were created, which represented the infected hosts and infected parasitoids. These two species were meant to represent the growth stage during which the host can no longer be considered “infectable” but is not yet a parasitoid or hyperparasitoid. The results were found to be very similar to those of the discrete models and the models are available upon request from


Parasitoids are a significant species when considering biological controls. Biological control refers to pest control of a certain species. Basically, the pest's natural enemy is introduced to keep the species population in control. A biological control must be host specific, meaning it has a small range of hosts it can infect. For example, a species many are somewhat familiar with is the Cane Toad, a biological control introduced in Australia to control pests which ended up eating anything and everything because it was not host specific. One other factor that plays into a hyperparasitoid-parasitoid-host system is host specificity. Host specificity is the collection of hosts a parasitoid can or will infect; the host range, in other words. Hyperparasitoids, which have only evolved in three orders of insecta; Hymenaptera (micro-wasps), Diptera (flies), and Coleoptera (beetles); are not very host-specific. Many times, they can infect an uninfected host, or if that host is already infected by a parasitoid, the hyperparasitoid can parasitize that primary parasitoid. As you may guess, hyperparasitoids can really throw off a biological control system. So that is one area in which hyperparasitoids play a big ecological role relevant to people.


Throughout this study, we progressed through a series of different predator-prey models. Predator-prey models originated from the Lotka-Volterra model of 1910 to determine interactions in continuous time. The Nicholson-Bailey model followed the Lotka-Volterra model, but changed the model to discretized time because it was more realistic and easier to measure. However, the model is unstable. Therefore, Hassell and May modified the Nicholson-Bailey model, by redefining the function in the system of difference equations, so that the equilibrium points were stable. Next, Hassell and May included a non-random search extension in which predators (or parasitoids) were modeled with a more specific searching pattern. We replicated both the random and non-random Hassel and May search models, and simulated both in continuous time. We also extended their model to include hyperparasitoids, creating a three species system. When we simulated the random search, the results were unstable. But when we simulated the non-random search pattern, the solutions were again stable.


  • Hassell, Michael P. The Dynamics of Arthropod Predator-prey Systems. Princeton, NJ: Princeton UP, 1978. Print.
  • Hassell, M. P. and R. M. May. "Stability in Insect Host-Parasite Models." Journal of Animal Ecology 42 (1973): 693 - 726.
  • Roberts, Larry S; Schmidt, Gerald D. Foundations of Parasitology. Edition 8. New York, New York: McGraw-Hill, 2009.


Here are the matlab scripts used in the creation of the plots for this project:

In modeling the random search for parasitoids:


In modeling the non-random search for parasitoids:


In modeling the random search for hyperparasitoids:


In modeling the non-random search for hyperparasitoids:


To test the speed of convergence for based off various variables


Modified version of the non-random search for hyperparasitoids to be used with the stability tester: