A Mathematical Investigation of Populations and Predator-Prey Dynamics

    • Introduction

In the late 1980s, research scientists working for the National Forest Service in the Sierra Nevada con-ducted a study to understand why a native population of mule deer had dramatically decreased in size from a 1950s census. The investigators believed that the initial reason for the decrease in the deer population was over-population of the deer habitat. However, the population continued to decrease even after it had been reduced to a size that the environment could support. The researchers wanted to nd out why the deer population never recovered.

Over the course of the study, the investigators noticed that the native mountain lion population had steadily increased while the deer population decreased. Mountain lions are natural predators of deer. It is generally believed that a predator species such as a mountain lion bene ts its prey population by keeping the prey population in balance with the environment. It does this by removing weak and old members of the prey population that would otherwise consume resources needed for younger, healthier members of the population to survive. In this way, the mountain lions (predators) are dependant upon the deer (prey) for their survival. The interactions of these species in uence the evolution of each population over time.

By studying the interactions of the deer and mountain lions, the researchers were able to deduce reasons why the deer population remained depressed years after its initial decline.

    • Modeling Individual Populations: the Logistic Equations

Suppose that these scientists wanted to study the populations of mountain lions and deer mathematically. One way to approach this task is to individually model each population to see how it changes with time (under a number of simplifying assumptions).

Assume the mountain lions are protected from hunting and have no natural predators. The mountain lions depend on the deer for food, but the deer population is nite. Thus, the amount of food available to the mountain lions is limited. Since the two populations are being studied separately, the mountain lion model must account for this constraint on its size without directly including the interactions of the mountain lions and the deer. A reasonable way to represent the change in the mountain lion population is with the logistic equation
dx
= r
1

x
x
(1)









dt


L


In this equation, x(t) is the size of the population (in dozens of animals) at any given time t (in days). The parameter r > 0 is the initial growth rate of the population and L > 0 is the carrying capacity of the population1. An underlying assumption of this model is that the population will grow exponentially in the absence of any external constraints on its size. We see this exponential behavior if we only include the rst


1Keep in mind that the accuracy of any prediction based on the logistic model depends upon whether the parameters r and L are constant.

1
term, rx(t), of the equation. Of course, a population is naturally limited in size by the available resources such as food or shelter. The second term in the equation, rx(t)2=L, is a correction factor that mimics the e ect of environmental constraints naturally limiting the population’s growth.

To model the deer population, it is reasonable to start from a logistic-type model because the deer are similarly limited by the resources in their environment. However, mountain lions are a natural predator of deer and are therefore a secondary means of restraining the deer population. It is sensible that the deer population model should include a predation term. A simple way to include the e ect of predation on the deer population is to modify the logistic equation to include a \harvesting" term H(x) that represents the number of deer killed by mountain lions. The modi ed logistic equation with harvesting becomes

dx
= r
1
x
x  H(x)
(2)








dt


L


where the meaning of x(t), r, and L are the same as in (1). A reasonable harvesting function could be

px2

H(x) = q + x2
(3)


The parameters p and q represent how skilled the mountain lions are at catching deer.

Section 2: Questions

    1. What are the units of the parameters r and L?

    2. Find the equilibrium solutions of (1) (i.e. the equation without harvesting). Then use separation of variables to nd the exact (non-equilibrium) solution to (1) (your answer should be reported in terms of the unspeci ed variables and parameters x, t, r, L, and some unknown initial population x(0) = x0).

    3. Suppose that the mountain lion population has a natural growth rate of r = 0:65 and a carrying
capacity of LM = 5:4.

        (a) If there are initially 6 mountain lions in the population2, use Euler’s method over the interval t 2 [0; 25] with step sizes h1 = 0:5, h2 = 0:1, and h3 = 0:01 to numerically solve the logistic equation (1). Plot the numerical solutions against the exact solution found in Question 2 above. (The plot should contain three numerical solution curves, the exact solution curve, and a legend). Discuss the plot. Speci cally address how the value of the step size h in uences the accuracy of the numerical solution.

        (b) The absolute error of a numerical approximation is de ned by

Abs. Error = jExact Solution    Approximate Solutionj

Use semilogy to plot the absolute errors of all numerical solutions found in part (a) (plot these error curves together on one plot. Be sure to include a legend). Describe the plot. In particular, speculate about why the error curves contain downward \spikes" around time t = 7.

    (c) Whenever a problem is solved numerically, it is important to make sure that the calculation carefully balances numerical \accuracy" { the correctness of the approximation obtained { with numerical \e ciency," which is the amount of time needed to complete the calculation. For the calculations in part (a) above, which step size might give the best balance of numerical accuracy and e ciency? Justify your response (the result of part (b) may help).


    • Hint: if we measure the mountain lion population in \dozens of mountain lions" then how you should write down this initial condition?



2
    4. Classify the di erential equation (2) (i.e. the equation with harvesting). Is it linear? What is its degree? Is it autonomous? Describe the physical meaning of autonomy or non-autonomy for this equation.

    5. Explore the behavior of the harvesting function (3) (a plot might help). What happens to H(x) as x becomes very large? What if x is close to 0? Does this make sense physically? Explain.

    6. Suppose the deer population also has a growth rate of r = 0:65 but has a carrying capacity of
LD = 8:1 (also let p = 1:2 and q = 1).

        (a) Modify the provided code dirfield.m in Matlab so that it computes the direction eld of the logistic equation with harvesting (2).

        (b) If there are initially 24 deer in the population3, use Euler’s method with a step size of h = 0:1 to numerically solve the logistic equation with harvesting (2) over the interval t 2 [0; 75]. Pick one other initial condition and use Euler’s method with the same h to numerically solve (2) over the same t interval.

        (c) Plot the Euler solutions from (b) and the direction eld of (2) together against t (include a legend for the solution curves). Discuss the plot. In particular, describe how the behavior of the solution curves changes with the initial condition and speculate about why this might be. (Hint: think about the locations and stability of equilibrium solutions.)

    • Modeling Population Interactions

Other than the harvesting term used to describe the e ect of predation on the deer population in (2), neither of the logistic equations (1) and (2) directly accounts for the interactions of the two populations. To more formally study how the populations respond to the pressures of both internal competition and external predator-prey interactions, the models must be adjusted to include the interactions of predator and prey and then solved simultaneously as a system. There are a number of ways to account for the interdependence of the two species.

3.1    The Lotka-Volterra System

One of the simplest models of predator-prey interactions is the Lotka-Volterra system

dx1
=  ax1 + bx1x2

dt





dx2
= cx2    dx1x2
(4)
dt






The Lotka-Volterra system is a result of the Balance Law: the net rate of change of a population is equal to the rate-in of members (birth/immigration) minus the rate-out of members (death/emigration). In this model, the variable x1 represents the predator population and x2 represents the prey population. These variables both reside in the rst quadrant of the x1x2-plane, which is also called the population quadrant. The positive parameters a, b, c, and d can be interpreted as follows:


a


predator mortality rate






b

predator attack rate/ conversion e  ciency (food into o spring)






c


prey growth rate






d

prey mortality rate/ predator attack rate











    • Hint: if we measure the deer population in \dozens of deer" then how you should write down this initial condition?

3
The cross terms, bx1x2 and dx1x2, in this model represent the interactions of the two species.

Notice that the predator population is a ected positively and the prey population is a ected negatively by interactions. In other words, the abundance of food (prey) promotes the predator’s growth rate but the presence of predators diminishes the prey’s growth rate. A basic assumption of this model is that, in the absence of interaction, each population will obey an exponential growth model wherein the predator population x01 = ax1 decays exponentially to zero while the prey population x02 = bx2 grows exponentially and without bound.

Section 3.1: Questions

    1. Classify the Lotka-Volterra system (4). Is it linear? Autonomous? What is its order?

    2. Analytically nd the nullclines and all equilibrium solutions of (4). DO NOT use speci c values for any of the parameters a, b, c, or d.

    3. Assign the parameters in (4) the values a = 1:5, b = 1:1, c = 2:5, and d = 1:4. Modify the provided script flow.m in Matlab so it computes the phase portrait of the system (4). Plot the following together in the x1x2-plane and label the curves appropriately (note: the \x1x2-plane" means that the x1 species is plotted on the x-axis and the x2 species is plotted on the y-axis):

Use flow.m to plot the phase portrait of (4) on the domain    1 < x1 < 6,    1 < x2 < 6.
Plot the nullclines of (4) on the domain    1 < x1 < 6,    1 < x2 < 6.

            ▪ Plot all nullclines that come from the x01 equation in red

            ▪ Plot all nullclines that come from the x02 equation in black Mark all equilibrium solutions with green.

    4. Classify the stability of the equilibrium solution (x1; x2) = (0; 0). What appears to be happening around the other equilibrium solution?

    5. Use ode45 in Matlab to simulate solutions to (4) starting from the initial condition x1(0) = 0:5,
x2(0) = 1:0 over the time interval t 2 [0; 20], then plot the following together:

Use flow.m to plot the phase portrait of (4) on the domain    1 < x1 < 6,    1 < x2 < 6.
Plot the phase-plane4 solution (x1; x2) on the domain    1 < x1 < 6,    1 < x2 < 6.

Does the solution curve behave as expected? Why or why not?

    6. Use ode45 in Matlab to simulate solutions to (4) starting from the initial condition x1(0) = 0:5,
x2(0) = 1:0 over the time interval t 2 [0; 20]. Plot the component curves5 x1(t) and x2(t) together against t (include a legend). Discuss these plots. Are the curves in phase or out of phase? What does this mean physically in terms of predator-prey interactions?









    • Note: the \phase-plane" solution is given by all sets of parametric points (x1(t); x2(t)). Generate the desired plot by plotting x2 against x1 in the x1x2-plane.

5Note: the \component curve" solutions are given by all pairs of points (t; x1(t)) and (t; x2(t)). Generate the desired plot by plotting x1 and x2 against t

4
3.2    The Logistic Predator-Prey

Recall that an underlying assumption of (4) is that both species will exhibit exponential behavior if there are no inter-species interactions (i.e. if the interaction terms x1x2 are excluded from the model). This assumption ignores the natural limits imposed on a prey population by its environment, such as nite food. A slightly more complicated but realistic way to model predator-prey interactions is to adjust the prey equation x02 to include these constraints, which gives rise to the Logistic Predator-Prey system
dx1
=  ax1 + bx1x2

dt





dx2
= c(1  kx2)x2    dx1x2
(5)
dt






Now the underlying assumption of this model is, in the absence of predation, the prey population will obey a logistic growth model instead of an exponential growth model. This means that

If the prey population is small, the rate of growth is proportional to its size

If the prey population is too large to be supported by its environment, the population will decrease As in (4), if the species do not interact the predator population will exhibit exponentially decaying solutions.

Section 3.2: Questions

    1. Analytically  nd the nullclines and equilibrium solutions of the Logistic Predator-Prey system (5).

    2. Assign the parameters in (5) the values a = 1:5, b = 1:1, c = 2:5, d = 1:4, and k = 0:5. Modify the provided script flow.m in Matlab so it computes the phase portrait of the system (5). Use ode45 in Matlab to simulate solutions to (5) starting from the initial condition x1(0) = 0:5, x2(0) = 1:0 over the time interval t 2 [0; 50] and then plot the following together:
Use flow.m to plot the phase portrait of (5) on the domain    1 < x1 < 6,    1 < x2 < 6.
Plot the nullclines and the equilibrium solutions of (5) on the domain    1 < x1 < 6,    1 < x2 < 6.
Plot the phase-plane solution (x1; x2) on the domain    1 < x1 < 6,    1 < x2 < 6.

Discuss the plot. What can you say about the stability of the equilibrium solutions? How does this in uence the solution curve?

    3. Assign the parameters in the Logistic Predator-Prey system (5) the values a = 1:5, b = 1:1, c = 2:5,

d = 1:4, and k = 0:5.  Use ode45 in Matlab to simulate solutions to (5) starting from the initial
condition x1(0) = 0:5, x2(0) = 1:0 over the time interval t 2 [0; 50]. Plot the component curves x1(t) and x2(t) together against t (include a legend). Discuss the solution curves. Are they periodic? Is there asymptotic behavior?

3.3    Model Comparison

Compare and contrast the Lotka Volterra (4) and the Logistic Predator Prey (5) models. What are some strengths and weaknesses of each model? Propose a modi cation to one or both of these models that might increase the accuracy of their predictions. Why/how does this improve the model?






5
    • References

Predators and PreyA Case of Imbalance Mountain Lions and the North Kings Deer Herd. Johnston Ridge Observatory | US Forest Service, www.fs.fed.us/psw/publications/Popular/mtnlions.html.

Mitotic Index and Cell Division, www.tiem.utk.edu/ gross/bioed/bealsmodules/predator-prey.html.

    • Project Guidelines

Your group will submit your project write-up on Canvas to the appropriate \Project Assignment" (you can nd these under the \assignments" tab in Canvas). Adhere to the following guidelines:

Do not put o nding a group (you must works in groups of 2-3). You should have a group set-up within one week of the project assignment due date.

Submit your project in a pdf format and submit ALL code used for your project (.nb les for Mathematica, .m les for Matlab, .py or ipynb for python). Code in Matlab, Mathematica, Maple, R, C, C++, Python or Julia is acceptable (Matlab is recommended). Code in Microsoft Word or Excel (or any other spreadsheet program) is not acceptable. All other languages need instructor permission (please ask as soon as possible).

Code may be included in the appendix if you wish. DO NOT submit anything on Canvas as a .zip le. Contents of .zip les will not necessarily be graded.

Have only ONE group member submit the project. Having multiple people in your group submit the project to Canvas will result in multiple grades, and we will take the LOWEST one.

Include the names and recitation section numbers of all group members working on the project on the cover page of the report.

When you submit the report to Canvas, please include each group member’s information (name, student number, and section number) in the comments. This allows us to quickly search for a student’s report.

Your report needs to accurately and consistently describe the steps you took to answer the posed questions.

This report should have the look and feel of a technical paper. Presentation and clarity are very important.

Here are some important items to remember:

Remember: you are to submit a complete report for this project. Documents submitted with numbered responses will be severely penalized.

Labs must be typed, including all equations (part of your learning experience is to learn how to use an equation editor). An exception can be made for lengthy calculations in the appendix, which may be hand written (as long as they are neat and clear), and minor labels on plots, arrows in the text, and a few subscripts.

Write your report in an organized and logical fashion. Section headers such as Introduction, Back-ground, Problem Statement, Calculations, Results, Conclusion, : : : are not mandatory but are highly recommended. They not only help you write your report, but help the reader navigate your paper.

Start with an introduction that describes what you will discuss in the body of the report. A brief summary of important concepts used in your discussion could be helpful here as well. Always introduce relevant equations that will be used or discussed in the report.



6
Always include units in your answers

DO include and label any plot that supports your conclusions or gives you insight in your investigations (these should be found in the body of your report). However, DO NOT use screen-shots of your  gures.

DO NOT include printouts of computer software screens or code snippets. You simply need to state which software you used in each step and what it did for you.

DO include the results of any calculations in the main body of your report. DO NOT include lengthy calculations in the main body of your paper (calculations should be included in a labeled appendix and should be referred to in the main body).

Your report does not have to be long. You need quality, not quantity, of work. Do not omit any important piece of information, but do not feel obligated to add any extras.

Summarize what you have accomplished in the conclusion. No new information or new results should appear in your conclusion. You should only review the highlights of what your wrote about in the body of the report.














































7