Predator Prey Graphing Assignment For Chemistry

The population density of collared lemmings (Dicrostonyx groenlandicus) is measured against the density of offspring produced by one of its predators, the Arctic fox (Vulpes lagopus). The data were collected in the High-Arctic tundra in Greenland after annual snowmelt.

Download the educator guide PDF, which includes questions to guide a class discussion on the graph characteristics and what the data show. Show the image on an overhead projector and lead the class in discussion or use the student handout before discussion.

Figure 1 from:

Gilg, Olivier et al., Cyclic dynamics in a simple vertebrate predator-prey community. Science. 2003, 302(5646), 866−868.

View article:

Science in the Classroom Annotated Paper:

Date Created 09/25/2017

Date Modified 10/12/2017

Predator-prey models are argubly the building blocks of the bio- and ecosystems as biomasses are grown out of their resource masses. Species compete, evolve and disperse simply for the purpose of seeking resources to sustain their struggle for their very existence. Depending on their specific settings of applications, they can take the forms of resource-consumer, plant-herbivore, parasite-host, tumor cells (virus)-immune system, susceptible-infectious interactions, etc. They deal with the general loss-win interactions and hence may have applications outside of ecosystems. When seemingly competitive interactions are carefully examined, they are often in fact some forms of predator-prey interaction in disguise.

A General Predator-Prey Model

Consider two populations whose sizes at a reference time \(t\) are denoted by \(x(t)\ ,\) \(y(t)\ ,\) respectively. The functions \(x\) and \(y\) might denote population numbers or concentrations (number per area) or some other scaled measure of the populations sizes, but are taken to be continuous functions. Changes in population size with time are described by the time derivatives \(\dot x \equiv dx/dt\) and \(\dot y \equiv dy/dt\ ,\) respectively, and a general model of interacting populations is written in terms of two autonomous differential equations \[ \dot x = x f(x,y) \] \[ \dot y = y g(x,y)\] (i.e., the time \(t\) does not appear explicitly in the functions \(x f(x,y)\) and \(y g(x,y)\)). The functions \(f\) and \(g\) denote the respective per capita growth rates of the two species. It is assumed that \( df(x, y)/dy<0 \) and \( dg(x, y)/dx>0. \) This general model is often called Kolmogorov's predator-prey model (Freedman 1980, Brauer and Castillo-Chavez 2000).

Lotka-Volterra Model

In 1926, the famous Italian mathematician Vito Volterra proposed a differential equation model to explain the observed increase in predator fish (and corresponding decrease in prey fish) in the Adriatic Sea during World War I. At the same time in the United States, the equations studied by Volterra were derived independently by Alfred Lotka (1925) to describe a hypothetical chemical reaction in which the chemical concentrations oscillate. The Lotka-Volterra model is the simplest model of predator-prey interactions. It is based on linear per capita growth rates, which are written as \[f= b-p y\] and \(g=r x-d\ .\)

  • The parameter \(b\) is the growth rate of species \(x\) (the prey) in the absence of interaction with species \(y\) (the predators). Prey numbers are diminished by these interactions: The per capita growth rate decreases (here linearly) with increasing \(y\ ,\) possibly becoming negative.
  • The parameter \(p\) measures the impact of predation on \(\dot x/x\ .\)
  • The parameter \(d\) is the death (or emigration) rate of species \(y\) in the absence of interaction with species \(x\ .\)
  • The term \(r x\) denotes the net rate of growth (or immigration) of the predator population in response to the size of the prey population.

The Prey-Predator model with linear per capita growth rates is \[\dot x = (b - p y) x\] (Prey) \[\dot y = (r x - d) y\] (Predators) This system is referred to as the Lotka-Volterra model: it represents one of the earliest models in mathematical ecology.

The system can be integrated directly. In particular, any solution \((x(t),y(t))\) of the system satisfies the identity \[ C = b \ln y(t) - p y(t) - r x(t) + d \ln x(t) \] for all \(t\ ,\) where the constant \(C= b \ln y(0) - p y(0) - r x(0) + d \ln x(0)\) is determined by initial conditions and system parameters. (Here \(\ln x\) denotes the natural logarithm of \(x\ ,\) etc.) The quantity on the right hand side of the identity above is referred to as a conservation law, since it is constant along any solution. Having a conserved quantity facilitates visualizing solutions. In the figure we draw level-set contours of the surface \[z = b \ln y - p y - r x + d \ln x\] in the first quadrant of the \(xy\)-plane. The contours describe solutions of the system determined by their initial data, and since they are closed curves, the solutions are periodic oscillations.

If \(b>0\ ,\) there are two equilibria, \(x=0, y=0\) (extinction), and \(x=d/r, y=b/p\) (co-existence), and the surface \[z = b \ln y - p y - r x + d \ln x\] has a single peak at the latter equilibrium. The contour lines in the figure describe the classic prey-predator cycles observed in ecological systems.

The model above has been derived independently in the following fields:

  • epidemics (Kermak and McKendrick 1927, 1932, 1933) (\(b=0\))
    • \(x\) are susceptible individuals and
    • \(y\) are infective individuals,
  • ecology (Lotka 1925, Volterra 1926)
    • \(x\) are prey and
    • \(y\) are predators,
  • combustion theory (Semenov 1935)
    • \(x\) and \(y\) are chemical radicals formed during H2 O2 combustion,
  • economics (Galbraith 2006)
    • \(x\) is the populace and
    • \(y\) is a predatory institution,

and numerous other studies from diverse disciplines.

Alfred James Lotka (1880 - 1949, USA, chemist, demographer, ecologist and mathematician) was born in Lviv (Lemberg), at that time situated in Austria, now in Ukraine. He came to the United States in 1902 and wrote a number of theoretical articles on chemical oscillations during the early decades of the twentieth century, and authored a book on theoretical biology (1925). He then left (academic) science and spent the majority of his working life at an insurance company (Metropolitan Life). In that capacity he became president of the PAA (the Population Association of America).

Vito Volterra (May 3, 1860 - October 11, 1940, Italian, mathematician and physicist, best known for his contributions to mathematical biology) was born in Ancona, into a very poor family. Volterra showed early promise in mathematics before attending the University of Pisa. His work is summarised in his book Theory of Functionals and of Integral and Integro-Differential Equations (1930). After World War I, Volterra turned his attention to the application of his mathematical ideas to biology, principally reiterating and developing the work of Pierre François Verhulst. The most famous outcome of this period is the Lotka-Volterra model. In 1922, he joined the opposition to the Fascist regime of Benito Mussolini and in 1931 refused to take a mandatory oath of loyalty. He was compelled to resign his university post and membership of scientific academies, and, during the following years, he lived largely abroad, returning to Rome just before his death.

Kermack-McKendrick Model

There is herd immunity in predation and in epidemics. It is convenient to frame this in terms of epidemiology where now we refer to the prey as being susceptibles and the predators as being infectives. The infection dynamics is depicted by the graph \(x\to y\to \ ,\) indicating that susceptibles can become infectives and that infectives can be removed from the process (e.g., through death, quarantine, or inoculation). Consider a time interval that is short compared to reproduction of the susceptible population, i.e., let \(b=0\ .\) If the initial susceptible population is so large that \(x(0)> d/r\ ,\) then we see from the second equation in the predator-prey model that initially \(\dot y>0\ ,\) which indicates that the infectives will initially more than replace themselves by passing on the infection. However, if this condition is not satisfied, the infective population will decrease. The critical value \(R\equiv rx(0)/d=1\) is referred to as being an epidemic threshold or a tipping point for the process.

A short calculation shows that \(x(t)\) converges to a constant, say \(x(t)\to x^*,\) where \(x^*\) can be found by solving the equation \(C=r x^* - d \ln x^*\ ,\) as shown in the figure. Surprisingly, this number is always greater than zero, which shows that some susceptibles will always survive!

This phenomenon, which is referred to as herd immunity is observed in practice; in fact, the number \(R\) is published regularly for various diseases and locales as an epidemic control measure. It reflects the fact that the susceptible population can be reduced to a level below which infectives will not increase. The model in this case is referred to as being the Kermack-McKendrick model of susceptible-infective interactions in epidemiology.

Jacob-Monod Model

Another approach to modeling the interaction between prey and predators was developed to account as well for organisms (such as bacteria) taking up nutrients. There is a limited uptake rate that such organisms are capable of, and the next model accounts for limited uptake rates. Suppose now that \(x\) denotes the size of a population of feeders and that they are feeding on a chemical species of concentration \(y\ .\) These are usually taken to represent concentrations of feeders and nutrients in solution rather than "head counts". The Jacob-Monod model is \[\dot x = \frac{V y}{K+y}x\] \[\dot y = -\frac 1Y \frac{V y}{K+y}x\] where

  • \(V\) is the uptake velocity,
  • \(K\) is the saturation constant, and
  • \(Y\) is the yield of \(x\) per unit \(y\) taken up.

Note that when \(y=K\ ,\) the uptake velocity is \(V/2\ ,\) half the maximum; in practice, \(y=K\) is taken as a tipping point: If \(y<K\ ,\) then uptake is ignored. This is the canonical model of nutrient uptake (nutrient \(y\) is taken up by species \(x\)), and it underlies many calculations in biology, microbiology, and food engineering (digestion, beer, etc.). This model was discovered independently in several diverse applications. It is akin to the Haldane-Briggs model and Michaelis-Menten model in biochemistry, the Jacob-Monod model in microbial ecology, and the Beverton-Holt model in fisheries. It serves as one of the important building blocks in studies of complex biochemical reactions and in ecology (Smith and Waltman 1997).

The quantity \(C=x +Y y\) is conserved (since the derivative with respect to time of the right hand side is zero) where \(C=x(0)+Yy(0)\ .\) Substituting this into the first equation gives \[\dot x =\frac{V (C-x)}{Y K+(C-x)}x\ .\] The solutions can be found by quadratures, and these show that if \(x(0)>0\ ,\) then \(x(t)\to C\) as \(t\to\infty\ ,\) at which point the nutrient has been depleted.

A typical use of this model is to describe a continuous-flow growth device, such as a chemostat, where there is continuous removal of nutrient and feeders and a continuous supply of fresh nutrient. The Jacob-Monod model is used to describe such a bacterial growth device, for example to determine conditions for a sustained dynamic equilibrium to exist by balancing growth due to uptake of nutrient with wash out of feeders.

Some predator-prey models use terms similar to those appearing in the Jacob-Monod model to describe the rate at which predators consume prey. More generally, any of the data in the Lotka-Volterra model can be taken to depend on prey density as appropriate for the system being studied. This is referred to as a functional reponse, an idea that is introduced and discussed by C. S. Holling (1959). Several different forms of functional response have been used in population models, but the Jacob-Monod form, also called the Holling type 2 form by ecologists, is one of the more common ones. Many other investigations of predator-prey models have involved functional responses. For example, \(r\) is replaced by \(r \ln(K/x)\) (Gompertz, 1825), \(r(K-x)/(K+\epsilon x),\) (Smith, 1963), and \(r\left((K/x)^g-1\right)\) with \(0< g \le 1\ ,\) (Rosenzweig, 1971). These and other functional responses are also discussed in May (1974). Such mechanisms in the Lotka-Volterra model can stabilize or destabilize the system, for example resulting in predator extinction or in co-existence of prey and predators. This is in contrast to the plurality of cycles predicted by the original Lotka-Volterra model.

Logistic Equation

An interesting case for \[\dot x =\frac{V (C-x)}{Y K+(C-x)}x\] is when \(V\) and \(YK\) are very large compared to the other data in the model, but with their ratio being of moderate size, say \(V/(YK)\approx r\ .\) Then we can ignore the second term in the denominator and get \[\dot x = r(C-x)x\ .\] This is referred to as the logistic equation. It, too, has arisen in various disciplines, but one of its first appearance was due to Verhulst in the mid 19th century who used it to correct certain deviations of Malthus's model (\(\dot x = r x\)) from certain human population data. The number \(C\) is now referred to as being the carrying capacity for the population; this corresponds to there being no remaining nutrient in the Jacob-Monod model.

The logistic equation can be solved in closed form by quadratures. This shows that \(x(t)\to C\) as \(t\to\infty\ ,\) if \(x(0)> 0\ .\)

Predation with Time Delays: Chaos in Ricker's Reproduction Equation

Time delays occur in biological systems, and they can produce complicated dynamics. To model age structure (and other time delays) in a system, we take the approach that was introduced by Euler in the 18th century. Let \(x(t)\) denote the number born at time \(t\ .\) Then \(x(t-a)\) denotes the number who were born \(a\) units ago. Suppose that there is a nonlinear effect such that the number of newborns \(a\) units ago who can participate in reproduction at time \(t\) is \(h(x(t-a))\) where the function \(h\) is called the reproduction function. This nonlinearity might be due to predation or environmental factors, as discussed earlier. Suppose those of age \(a\) have fertility \(m(a)\ .\) Then the equation below shows how many will be born at time \(t\ .\) It is obtained by adding up across all ages those who can participate weighted by their fertility: \[ x(t) = x_0(t) + \int_0^t m(a) h(x(t-a))\, da\ ,\] where \(x_0(t)\) denotes the lingering contributions to later births of the initial population. In the case where \(h(x)= x\) is a linear function, this equation is referred to as being the renewal equation, which is widely used in demography. Its solution can be found using Laplace transforms (Keyfitz and Flieger 1971). However, if \(h\) is a nonlinear function, then that method of solution is not available.

The fertility function \(m(a)\) describes the fertility of those of age \(a\ .\) If \(m\) is constant, then this equation is equivalent to a differential equation. The other extreme occurs when fertility is focused all at one age, as for salmon or cicadas. In this case, suppose that all fertility is due to those of age \(a^*\ .\) We write \(m(a)=r\,\delta(a-a^*)\) where \(\delta\) is Dirac's delta function and \(r\) is the reproduction rate. The equation becomes \[x(t) = x_0(t)+r h(x(t-a^*))\ .\] Let us focus only on the births at the times \(a^*, 2a^*, 3a^*,\ldots\ ;\) we define \(x_n=x(na^*)\) with \(x_1=x_0(a^*)\ .\) Substituting this in the equation above gives the recursion \[x_{n+1}= r h(x_n)\ ,\] for \(n=1,2,3,4,\ldots\ .\) In the case of Malthus's model, \(h(x) = x\ ,\) and the solution is a simple geometric progression \(x_n = r^n x_0\ .\) However, if \(h\) is as in Verhulst's case, where \(h(x) = x ( C - x)_+ \equiv \max(C-x,0),\) the recursion becomes \[x_{n+1}=r x_n(C-x_n)_+\ ,\] whose iterates can exhibit quite wild behavior. This was first noticed by the ecologist W.E. Ricker in the 1950's who used the function \(h(x)=x\exp{(-x)}\) in studies of the dynamics of fisheries, although his work was largely ignored at the time. This reproduction function accounts for cannibalism (self predation) in that if the population is small, the model looks like Malthus's, but for large populations, reproduction is strongly suppressed. The work was later rediscovered by Robert May, who stimulated the now prominent area of chaos (Gleick 1987).

Chaotic dynamics can be illustrated by a simple computer simulation using Ricker's model. Since the maximum of \(h\) is 1, we need only consider values of \(x\in [0, r]\ .\) We divide this interval up into 300 equal bins. The parameter \(r\) is fixed at one of 300 equally spaced numbers in \(r\in [0,20]\ .\) The recursion is \[x_{n+1}=r x_ne^{-x_n}\ .\] Beginning with \(x_1=1\ ,\) we determine the sequence \(\{x_{50},x_{51},x_{52},x_{53},\dots,x_{300}\}\) and record in an array, say \(D(r,x)\ ,\) the number of sequence elements that are in each bin. We only start counting after 50 initial iterates to avoid transients. The result is a surface \(D(r,x)\) that describes the dynamics. We plot this record looking from above, as shown in the figure (Hoppensteadt and Hyman 1977). The solution always settles into some structure, which might be highly complex. At first, there is a unique stable state. Near \(r=8\) this state bifurcates into an oscillation of period 2. Near \(r=12\) the period 2 solution bifurcates into a period 4 solution, etc. By \(r=18\) the solution set is quite complicated. This is in the range of chaotic dynamics. When the fertility function \(m\) has broader support, the solution set is more complicated still.


  • F. Brauer and C. Castillo-Chavez, Mathematical Models in Population Biology and Epidemiology, Springer-Verlag, Heidelberg, 2000.
  • H. I. Freedman, Deterministic Mathematical Models in Population Ecology. New York: Marcel Dekker, 1980.
  • J.K. Galbraith, The Predator State, Mother Jones, May/June 2006.
  • J. Gleick, Chaos: The Making of a New Science, Viking Press, New York, 1987.
  • C. S. Holling, The characteristics of simple type of predation and parasitism, Canadian Entomologist 91 (1959), 385-398.
  • F.C. Hoppensteadt, J.M. Hyman, Periodic solutions of a logistics difference equation, SIAM J. Appl. Math. 58(1977)73-81.
  • O. Kermack, A.G. McKendrick, Proc. Roy. Soc. A, 115(1927)700-721, 138(1932)55-83, 141(1933)94-122.
  • N. Keyfitz, W. Flieger, Populations: Fact and Methods of Demography, W.H.Freeman, San Francisco, 1971.
  • M. Kot, Elements of Mathematical Ecology, Cambridge University Press, 2001
  • A.J. Lotka, Elements of physical biology. Williams and Wilkins, Baltimore, 1925.
  • R. M. May, Stability and Complexity in Model Ecosystems, Princeton U. Press, NJ, 1974.
  • N.N. Semenov, Chemical Kinetics and Chain Reactions, Clarendon, Oxford, 1935.
  • H. Smith, P. Waltman, The Mathematical Theory of Chemostats, Cambridge U. Press, 1997.
  • V. Volterra, Variazioni e fluttuazioni del numero d'individui in specie animali conviventi. Mem. R. Accad. Naz. dei Lincei 2(1926)31�113

Internal references

  • Jeff Moehlis, Kresimir Josic, Eric T. Shea-Brown (2006) Periodic orbit. Scholarpedia, 1(7):1358.
  • Philip Holmes and Eric T. Shea-Brown (2006) Stability. Scholarpedia, 1(10):1838.

External Links

See also

Dynamical Systems, Periodic Orbit, Population Dynamics, Phase Space, Stability

Figure 1: Periodic activity generated by the Predator-Prey model.

Figure 2: Prey-Predator dynamics as described by the level curves of a conserved quantity. The arrows describe the velocity and direction of solutions. In this simulation, the data are \(d=r=b=d=1\ .\) There are equilibria at \(x=1, y=1\) and at \(x=0, y=0\ .\)

Figure 3: Kermack-McKendrick model of propagation of infectious disease. The tipping point is at \(x = 1.0\ .\) If \(x(0)\) is above this value, an epidemic will ensue. The severity can be estimated by tracing the curve emanating from \(x(0)\) until it converges to the horizontal axis. This indicates the size of the susceptible population that avoids the infection. Note that the horizontal axis as drawn here begins at x = 0.1 to avoid the logarithmic singularity.

Figure 4: Dynamics of Ricker's population.

0 thoughts on “Predator Prey Graphing Assignment For Chemistry

Leave a Reply

Your email address will not be published. Required fields are marked *