
Extensions to The Dynamics of Arthropod PredatorPrey SystemsFrom MathBioReview by Justin Baacke, Nicole Look and Brit Schneiders, May 2012. The following wiki reproduces the results presented by one of the insect hostparasite models developed in Arthropod PredatorPrey Systems by M. P. Hassell. More specifically, the article below first demonstrates the analytical and computational analysis of the model for nonrandom 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 nonrandom search, and also incorporating hyperparasitoidism. The work shown below can be attributed to Hassell unless otherwise specified. ContentsExecutive SummaryThe nonrandom search model presented below, as well as other models that describe species interactions, developed from modifying the LotkaVolterra predatorprey model of 1910. Although the LotkaVolterra model simplistically describes the interactions between two species, it only considers continuous time. Nicholson and Bailey (1935) resolved this problem by creating a predatorprey 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 NonRandom search. Random search is an easier model to create, but nonrandom search is a logical modification because it is seen all around in nature. In the following article, we reproduce both the random and nonrandom search models in discrete time. Then we also use an agentbased model to simulate what might happen for both random and nonrandom 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 twospecies system in a threespecies 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 twospecies and threespecies systems. Then, when we implemented nonrandom search, we found stable results in discretized and continuous time, for both parasitoids vs. hosts, and hyperparasitoids vs. parasitoids vs. hosts. Biological ContextMathematically, parasitoids are a good, simple model of predators when describing predatorprey 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 onetoone 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: 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 freeliving 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):
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:
Another biological concept that we will explore is that of hyperparasitoids. Hyperparasitoids, in terms of a predatorprey 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 wellstudied at the time Hassell publicized his book on Arthropod PredatorPrey Systems. HistoryPredatorprey interactions were originally explained by the LotkaVolterra 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 LotkaVolterra model to incorporate discrete time. More specifically, the model describes the population dynamics of a coupled hostparasite 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 hostparasite interactions. In addition, it yields diverging oscillations after perturbations. Most species populations do not tend toward extinction. For a more indepth explanation of the NicholsonBailey model, please consult “The Dynamics of Arthropod PredatorPrey Systems.” To add stability and introduce nonrandom search, Hassell and May (1973) then modified the NicholsonBailey model. By considering a nonhomogeneous distribution of prey and a nonrandom 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 ResultsThe majority of insect hostparasite models follow a general form: (1) where is the number of host survivors once parasites have searched for hosts that results in a subsequent generation of parasite progeny. Modifications of the general model can be made by changing the function , 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 hostparasite interactions is affected by nonrandom search. This proposition contradicted the previous hostparasite 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 nonrandom search pattern and take advantage of the highdensity 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 PredatorPrey 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) (2) with the newly defined variables: Fraction of host population in patch Fraction of parasite population in patch Number of patches Prey rate of increase Equation (1) can then be modified as: (3) The summation of the considered space allows for prey to nonuniformly occupy any of the patches, such that parasites and hosts can be distributed into n areas specified by and . 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 and function (2) simplifies to the NicholsonBailey model.
(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 and : (5) provided that . 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 , the perturbed populations can be written as (6) where and are small initial perturbations. Given the total number of prey attacked by predators, the dynamics of the perturbations can be studied by Taylor expansion of (7) about the equilibrium. After discarding higher order terms (8) where and and the partial derivatives are evaluated at the equilibrium point. Since the fraction of unparasitized hosts most likely decreases when the number of parasites increases and increases when the number of hosts increases, it is expected that Standard methods for analyzing stability transform the solution of the linear difference equations into the form of (9) By substitution of equations (9) into equations (8), the expanded equations can be rewritten as (10) Equations (10) are only consistent when leading to the characteristic equation (11) which can be solved as (12). The overall stability criterion is :
Applying the stability criterion to solution (12) leads to the overall stability criterion of , such that and (13). It can be noted that
More specific for the nonrandom search model, where . In this case, the equilibrium populations and found in equation (5) can be applied to equation (8) to yield and as (14). The overall stability criterion (13) leads to the exact condition for stability (15). It is important to note that stability is most affected by:
Furthermore, the scenario analyzed by Hassell and May considered one patch with a high density fraction, , of the population and the remaining "n1" patches containing as a density fraction of the prey. This allows for the introduction of an aggregation index which is defined as (16) where c is a normalizing constant. Here is considered a measure of the distribution of predators. When a random search pattern exists, and when , 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:
Nonrandom 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 ModelComputational 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 multipatch 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. 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: 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. 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 nonconvergence from .6 to .85 and then we resume the relatively quick convergence which is likely due to the system dying off. 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. 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. 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. Agent Based Models In order to extend our model into continuous space we also analyzed the effect of a random search path, versus a nonrandom 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 Justin.baacke@colorado.edu Extension: HyperparasitoidsHyperparasitoids, 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 predatorprey system, or in this case, a parasitoidhost 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 predatorprey 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! 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. 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. 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 nonrandom 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 Justin.baacke@colorado.edu ApplicationsParasitoids 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 hyperparasitoidparasitoidhost 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 (microwasps), Diptera (flies), and Coleoptera (beetles); are not very hostspecific. 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. ConclusionThroughout this study, we progressed through a series of different predatorprey models. Predatorprey models originated from the LotkaVolterra model of 1910 to determine interactions in continuous time. The NicholsonBailey model followed the LotkaVolterra 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 NicholsonBailey model, by redefining the function in the system of difference equations, so that the equilibrium points were stable. Next, Hassell and May included a nonrandom search extension in which predators (or parasitoids) were modeled with a more specific searching pattern. We replicated both the random and nonrandom 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 nonrandom search pattern, the solutions were again stable. References
CodeHere are the matlab scripts used in the creation of the plots for this project: In modeling the random search for parasitoids: In modeling the nonrandom search for parasitoids: In modeling the random search for hyperparasitoids: In modeling the nonrandom search for hyperparasitoids: To test the speed of convergence for based off various variables Modified version of the nonrandom search for hyperparasitoids to be used with the stability tester: 