predator_prey_lab.sagews

Modeling population interactions in ecosystems

This lab's theme is based on a true story that illustrates the complex interdependence of species in many ecosystems. In the late 19th century the U.S. citrus industry was nearly destroyed by the accidental introduction of an insect from Australia (the "scale-insect") that turned out to be a major pest. One of the many efforts the industry made to combat the pest was to introduce a species of beetle known to be a natural predator of the Australian scale-insect. Over time the industry found this strategy unsatisfactory and decided to spray DDT, a toxic pesticide, to get rid of both species. To their surprise, this made the problem worse. The population of the pest seemed to increase overall, while the population of the beetles decreased!

Source of story and modeling idea: Mathematical Modelling with Case Studies, by Belinda Barnes and Glenn Fulford.

Overall objectives:
In this lab we will build upon the concepts developed in our recent lab on rate of change models for different growth processes. This time we'll conjecture and design a model for the rate of change of two populations that interact with each other -- the Australian scale-insect (prey) and the beetle (predator).

We make the following two assumptions which, though inaccurate, lead to remarkably good real-world models:

1. The prey population is kept in check only by the predator. In their absence the prey population would grow exponentially to infinity.
2. The predator population tends to die out on its own. It only grows in the presence of the prey.

Let's use the following units and notation:

• $t$ = time in days
• $S(t)$ = scale-insect population (number per unit area)
• $B(t)$ = beetle population (number per unit area)

(I) Start with simple case - assume each population is on its own
(there is no interaction between them)

1. The assumptions directly tell us the approximate shape of the $S$ vs. $t$ and $B$ vs. $t$ graphs. Sketch graphs illustrating these qualitative behaviors. Explain or justify your graphs.
   Can you tell what the graph of graph $B$ vs. $S$ looks like?

• Suppose we have the following observational data showing the population change in one day for each species (in isolation from the other):
          The scale-insect goes from 500 to 540 per unit area.
          The beetles go from 50 to 49 per unit area.
  Write a mathematical statement of the rate of change of $S(t)$ and $B(t)$.
          Hint: Rate of change = derivative.
          $S^\prime(t)$ = negative term, causes $S$ to decrease.
          $S^\prime(t)$ = positive term, causes $S$ to increase.

  Some examples:

  $S^\prime=1$, OR $~S^\prime=t$, OR $~S^\prime=e^t$, OR $~S^\prime=S$, OR ...

  Clearly each of these corresponds to an increasing function $S(t)$.    [Why?]
  In that sense, each is a plausible model for the scale-insect population.
  But, which of them does the "best" job? Why?
  Write an analogous model for $B^\prime$.

(II) More realistic (& interesting) case - the predator and prey populations interact

Start with your simple model above and think about modifying it to account for the interaction between predator and prey. You want to come up with rate of change equations that look like:
$S^\prime(t)$ = stuff,   $B^\prime(t)$ = other-stuff.

Things to remember: $S^\prime(t)$ = negative term, causes $S(t)$ to decrease;
S$^\prime(t)$ = positive term, causes $S(t)$ to increase;
$S^\prime(t)$ = negative + positive term, causes $S(t)$ to respond to two opposing effects.

Here are some explorations to try:

1. Suppose your simple model for the scale-insects is $S^\prime = 0.08 S$.
   If you add a negative term to the right that depends on $B$, it will model the effect of beetles eating the scale-insects. Conjecture some possible forms for this 2nd term. Some ideas: $-B$, $~~-10B$, $~~-B^2$, $~-\sqrt{B}$, $~~-e^B$, $~~-S\cdot B$, $~~-0.5\;S\;B$, ...     What are your thoughts?
2. It turns out that a good model for the predator-prey effect is to have a term like $0.0002\;S\;B$ in your rate of change equations. Of course, I've put the "0.0002" there just to show you a specific example - it should be replaced by an appropriate numerical value for each specific application. With this in mind, write the general form of the $S^\prime(t)$ and $B^\prime(t)$ equations in our model - use symbols like a, b, c, etc., wherever you need constant coefficients:
       $S^\prime(t) =$ term1 + term2.
       $B^\prime(t) =$ term3 + term4.
   See page 240, exercise 35 in your textbook for a very similar model.
3. Now for the fun part! Simulate a predator-prey model using Sage. Let's use the following model:
       $S^\prime(t) = 0.08 S - 0.001 S B$
       $B^\prime(t) = -0.02 B + 0.00002 S B$

   We'll carry out the simulation for 3 different starting scenarios:
       (i) At t=0 there are 500 scale-insects and 50 beetles.
       (ii) At t=0 there are 500 scale-insects and 100 beetles.
       (iii) At t=0 there are 1000 scale-insects and 81 beetles.
   Consult with Anand for help on how to use Sage to solve this model.

   Discuss  what your simulations reveal. In particular, discuss what types of patterns you see in the predator and prey populations, and whether/ how they depend on the initial population of each species. For each simulation, use your graphical solutions to estimate an average value of the predator and prey populations. A very rough approximation of the averages is sufficient.

4. Effect of DDT: From a calculus point of view, the pesticide essentially contributes a new "death term" in the respective equations of each species. A question to ask yourself is, what is a reasonable form for this term: should it be constant, or should it depend on B or S? should the term be the same, or different, for each species?

   Run the grapher simulation again after adding the following new terms to your previous model:
       S'(t) = previous terms - 0.01 S.
       B'(t) = previous terms - 0.01 B.

   Comment on what your results reveal about the effect of using DDT. Address any changes you see in the population patterns, and in the average values.

# Systems is: S'= 0.08*S - 0.001*S*B; B'= -0.02*B + 0.00002*S*B; S(0)=500, B(0)=50

s, b, t = var('s b t')
de1 = 0.08*s - 0.001*s*b
de2 = -0.02*b + 0.00002*s*b
P = desolve_system_rk4 ([de1, de2], [s, b], ics=[0, 500, 50], ivar=t, end_points=[0,200] )
Q = [ [j,k] for i, j, k in P]
#list_plot(Q, axes_labels=['$S$', '$B$'])
line(Q, axes_labels=['$S$', '$B$'], color='green', thickness=2, title='Phase plot')

# The following creates plots of S vs t and B vs t:
QS = [ [i,j] for i, j, k in P]
QB = [ [i,k] for i, j, k in P]
SS = line(QS, axes_labels=['$t$ (days)', '$P(t)$'], legend_label='Scale insects', color='blue', thickness=2)
SB = line(QB, axes_labels=['$t$ (days)', '$P(t)$'], legend_label='Beetles', color='green', thickness=2)
show(SS+SB)