# 9.6: Predator-Prey Systems - Mathematics

Learning Objectives

• Describe the concept of environmental carrying capacity in the logistic model of population growth.
• Draw a direction field for a logistic equation and interpret the solution curves.
• Solve a logistic equation and interpret the results.

Differential equations can be used to represent the size of a population as it varies over time. We saw this in an earlier chapter in the section on exponential growth and decay, which is the simplest model. A more realistic model includes other factors that affect the growth of the population. In this section, we study the logistic differential equation and see how it applies to the study of population dynamics in the context of biology.

## Population Growth and Carrying Capacity

To model population growth using a differential equation, we first need to introduce some variables and relevant terms. The variable (t). will represent time. The units of time can be hours, days, weeks, months, or even years. Any given problem must specify the units used in that particular problem. The variable (P) will represent population. Since the population varies over time, it is understood to be a function of time. Therefore we use the notation (P(t)) for the population as a function of time. If (P(t)) is a differentiable function, then the first derivative (frac{dP}{dt}) represents the instantaneous rate of change of the population as a function of time.

In Exponential Growth and Decay, we studied the exponential growth and decay of populations and radioactive substances. An example of an exponential growth function is (P(t)=P_0e^{rt}.) In this function, (P(t)) represents the population at time (t,P_0) represents the initial population (population at time (t=0)), and the constant (r>0) is called the growth rate. Figure (PageIndex{1}) shows a graph of (P(t)=100e^{0.03t}). Here (P_0=100) and (r=0.03).

We can verify that the function (P(t)=P_0e^{rt}) satisfies the initial-value problem

[ dfrac{dP}{dt}=rP]

with (P(0)=P_0.)

This differential equation has an interesting interpretation. The left-hand side represents the rate at which the population increases (or decreases). The right-hand side is equal to a positive constant multiplied by the current population. Therefore the differential equation states that the rate at which the population increases is proportional to the population at that point in time. Furthermore, it states that the constant of proportionality never changes.

One problem with this function is its prediction that as time goes on, the population grows without bound. This is unrealistic in a real-world setting. Various factors limit the rate of growth of a particular population, including birth rate, death rate, food supply, predators, and so on. The growth constant (r) usually takes into consideration the birth and death rates but none of the other factors, and it can be interpreted as a net (birth minus death) percent growth rate per unit time. A natural question to ask is whether the population growth rate stays constant, or whether it changes over time. Biologists have found that in many biological systems, the population grows until a certain steady-state population is reached. This possibility is not taken into account with exponential growth. However, the concept of carrying capacity allows for the possibility that in a given area, only a certain number of a given organism or animal can thrive without running into resource issues.

Definition: Carrying Capacity

The carrying capacity of an organism in a given environment is defined to be the maximum population of that organism that the environment can sustain indefinitely.

We use the variable (K) to denote the carrying capacity. The growth rate is represented by the variable (r). Using these variables, we can define the logistic differential equation.

Definition: Logistic Differential Equation

Let (K) represent the carrying capacity for a particular organism in a given environment, and let (r) be a real number that represents the growth rate. The function (P(t)) represents the population of this organism as a function of time (t), and the constant (P_0) represents the initial population (population of the organism at time (t=0)). Then the logistic differential equation is

[dfrac{dP}{dt}=rPleft(1−dfrac{P}{K} ight). label{LogisticDiffEq}]

The logistic equation was first published by Pierre Verhulst in (1845). This differential equation can be coupled with the initial condition (P(0)=P_0) to form an initial-value problem for (P(t).)

Suppose that the initial population is small relative to the carrying capacity. Then (frac{P}{K}) is small, possibly close to zero. Thus, the quantity in parentheses on the right-hand side of Equation ef{LogisticDiffEq} is close to (1), and the right-hand side of this equation is close to (rP). If (r>0), then the population grows rapidly, resembling exponential growth.

However, as the population grows, the ratio (frac{P}{K}) also grows, because (K) is constant. If the population remains below the carrying capacity, then (frac{P}{K}) is less than (1), so (1−frac{P}{K}>0). Therefore the right-hand side of Equation ef{LogisticDiffEq} is still positive, but the quantity in parentheses gets smaller, and the growth rate decreases as a result. If (P=K) then the right-hand side is equal to zero, and the population does not change.

Now suppose that the population starts at a value higher than the carrying capacity. Then (frac{P}{K}>1,) and (1−frac{P}{K}<0). Then the right-hand side of Equation ef{LogisticDiffEq} is negative, and the population decreases. As long as (P>K), the population decreases. It never actually reaches K because (frac{dP}{dt}) will get smaller and smaller, but the population approaches the carrying capacity as (t) approaches infinity. This analysis can be represented visually by way of a phase line. A phase line describes the general behavior of a solution to an autonomous differential equation, depending on the initial condition. For the case of a carrying capacity in the logistic equation, the phase line is as shown in Figure (PageIndex{2}).

This phase line shows that when (P) is less than zero or greater than (K), the population decreases over time. When (P) is between (0) and (K), the population increases over time.

Example (PageIndex{1}): Examining the Carrying Capacity of a Deer Population

Let’s consider the population of white-tailed deer (Odocoileus virginianus) in the state of Kentucky. The Kentucky Department of Fish and Wildlife Resources (KDFWR) sets guidelines for hunting and fishing in the state. Before the hunting season of 2004, it estimated a population of 900,000 deer. Johnson notes: “A deer population that has plenty to eat and is not hunted by humans or other predators will double every three years.” (George Johnson, “The Problem of Exploding Deer Populations Has No Attractive Solutions,” January 12,2001, accessed April 9, 2015)

This observation corresponds to a rate of increase (r=dfrac{ln (2)}{3}=0.2311,) so the approximate growth rate is 23.11% per year. (This assumes that the population grows exponentially, which is reasonable––at least in the short term––with plentiful food supply and no predators.) The KDFWR also reports deer population densities for 32 counties in Kentucky, the average of which is approximately 27 deer per square mile. Suppose this is the deer density for the whole state (39,732 square miles). The carrying capacity (K) is 39,732 square miles times 27 deer per square mile, or 1,072,764 deer.

1. For this application, we have (P_0=900,000,K=1,072,764,) and (r=0.2311.) Substitute these values into Equation ef{LogisticDiffEq} and form the initial-value problem.
2. Solve the initial-value problem from part a.
3. According to this model, what will be the population in (3) years? Recall that the doubling time predicted by Johnson for the deer population was (3) years. How do these values compare?

Suppose the population managed to reach 1,200,000 What does the logistic equation predict will happen to the population in this scenario?

Solution

a. The initial value problem is

[ dfrac{dP}{dt}=0.2311P left(1−dfrac{P}{1,072,764} ight),,,P(0)=900,000. onumber]

b. The logistic equation is an autonomous differential equation, so we can use the method of separation of variables.

Step 1: Setting the right-hand side equal to zero gives (P=0) and (P=1,072,764.) This means that if the population starts at zero it will never change, and if it starts at the carrying capacity, it will never change.

Step 2: Rewrite the differential equation and multiply both sides by:

[ egin{align*} dfrac{dP}{dt} =0.2311Pleft(dfrac{1,072,764−P}{1,072,764} ight) [4pt] dP =0.2311Pleft(dfrac{1,072,764−P}{1,072,764} ight)dt [4pt] dfrac{dP}{P(1,072,764−P)} =dfrac{0.2311}{1,072,764}dt. end{align*}]

Step 3: Integrate both sides of the equation using partial fraction decomposition:

[ egin{align*} ∫dfrac{dP}{P(1,072,764−P)} =∫dfrac{0.2311}{1,072,764}dt [4pt] dfrac{1}{1,072,764}∫ left(dfrac{1}{P}+dfrac{1}{1,072,764−P} ight)dP =dfrac{0.2311t}{1,072,764}+C [4pt] dfrac{1}{1,072,764}left(ln |P|−ln |1,072,764−P| ight) =dfrac{0.2311t}{1,072,764}+C. end{align*} ]

Step 4: Multiply both sides by 1,072,764 and use the quotient rule for logarithms:

[ln left|dfrac{P}{1,072,764−P} ight|=0.2311t+C_1. onumber]

Here (C_1=1,072,764C.) Next exponentiate both sides and eliminate the absolute value:

[ egin{align*} e^{ln left|dfrac{P}{1,072,764−P} ight|} =e^{0.2311t + C_1} [4pt] left|dfrac{P}{1,072,764 - P} ight| =C_2e^{0.2311t} [4pt] dfrac{P}{1,072,764−P} =C_2e^{0.2311t}. end{align*}]

Here (C_2=e^{C_1}) but after eliminating the absolute value, it can be negative as well. Now solve for:

[ egin{align*} P =C_2e^{0.2311t}(1,072,764−P) [4pt] P =1,072,764C_2e^{0.2311t}−C_2Pe^{0.2311t} [4pt] P + C_2Pe^{0.2311t} = 1,072,764C_2e^{0.2311t} [4pt] P(1+C_2e^{0.2311t} =1,072,764C_2e^{0.2311t} [4pt] P(t) =dfrac{1,072,764C_2e^{0.2311t}}{1+C_2e^{0.23 onumber11t}}. end{align*}]

Step 5: To determine the value of (C_2), it is actually easier to go back a couple of steps to where (C_2) was defined. In particular, use the equation

[dfrac{P}{1,072,764−P}=C_2e^{0.2311t}. onumber]

The initial condition is (P(0)=900,000). Replace (P) with (900,000) and (t) with zero:

[ egin{align*} dfrac{P}{1,072,764−P} =C_2e^{0.2311t} [4pt] dfrac{900,000}{1,072,764−900,000} =C_2e^{0.2311(0)} [4pt] dfrac{900,000}{172,764} =C_2 [4pt] C_2 =dfrac{25,000}{4,799} [4pt] ≈5.209. end{align*}]

Therefore

[ egin{align*} P(t) =dfrac{1,072,764 left(dfrac{25000}{4799} ight)e^{0.2311t}}{1+(250004799)e^{0.2311t}}[4pt] =dfrac{1,072,764(25000)e^{0.2311t}}{4799+25000e^{0.2311t}.} end{align*}]

Dividing the numerator and denominator by 25,000 gives

[P(t)=dfrac{1,072,764e^{0.2311t}}{0.19196+e^{0.2311t}}. onumber]

Figure is a graph of this equation.

c. Using this model we can predict the population in 3 years.

[P(3)=dfrac{1,072,764e^{0.2311(3)}}{0.19196+e^{0.2311(3)}}≈978,830,deer onumber]

This is far short of twice the initial population of (900,000.) Remember that the doubling time is based on the assumption that the growth rate never changes, but the logistic model takes this possibility into account.

d. If the population reached 1,200,000 deer, then the new initial-value problem would be

[ dfrac{dP}{dt}=0.2311P left(1−dfrac{P}{1,072,764} ight), , P(0)=1,200,000. onumber]

The general solution to the differential equation would remain the same.

[ P(t)=dfrac{1,072,764C_2e^{0.2311t}}{1+C_2e^{0.2311t}} onumber]

[ dfrac{P}{1,072,764−P}=C_2e^{0.2311t}. onumber]

Substituting the values (t=0) and (P=1,200,000,) you get

[ egin{align*} C_2e^{0.2311(0)} =dfrac{1,200,000}{1,072,764−1,200,000} [4pt] C_2 =−dfrac{100,000}{10,603}≈−9.431.end{align*}]

Therefore

[ egin{align*} P(t) =dfrac{1,072,764C_2e^{0.2311t}}{1+C_2e^{0.2311t}} [4pt] =dfrac{1,072,764 left(−dfrac{100,000}{10,603} ight)e^{0.2311t}}{1+left(−dfrac{100,000}{10,603} ight)e^{0.2311t}} [4pt] =−dfrac{107,276,400,000e^{0.2311t}}{100,000e^{0.2311t}−10,603} [4pt] ≈dfrac{10,117,551e^{0.2311t}}{9.43129e^{0.2311t}−1} end{align*}]

This equation is graphed in Figure (PageIndex{5}).

## Solving the Logistic Differential Equation

The logistic differential equation is an autonomous differential equation, so we can use separation of variables to find the general solution, as we just did in Example (PageIndex{1}).

Step 1: Setting the right-hand side equal to zero leads to (P=0) and (P=K) as constant solutions. The first solution indicates that when there are no organisms present, the population will never grow. The second solution indicates that when the population starts at the carrying capacity, it will never change.

Step 2: Rewrite the differential equation in the form

[ dfrac{dP}{dt}=dfrac{rP(K−P)}{K}.]

Then multiply both sides by (dt) and divide both sides by (P(K−P).) This leads to

[ dfrac{dP}{P(K−P)}=dfrac{r}{K}dt.]

Multiply both sides of the equation by (K) and integrate:

[ ∫dfrac{K}{P(K−P)}dP=∫rdt. label{eq20a}]

The left-hand side of this equation can be integrated using partial fraction decomposition. We leave it to you to verify that

[ dfrac{K}{P(K−P)}=dfrac{1}{P}+dfrac{1}{K−P}.]

Then the Equation ef{eq20a} becomes

[ ∫dfrac{1}{P}+dfrac{1}{K−P}dP=∫rdt]

[ ln |P|−ln |K−P|=rt+C]

[ ln ∣dfrac{P}{K−P}∣=rt+C.]

Now exponentiate both sides of the equation to eliminate the natural logarithm:

[ e^{ln ∣dfrac{P}{K−P}∣}=e^{rt+C}]

[ ∣dfrac{P}{K−P}∣=e^Ce^{rt}.]

We define (C_1=e^c) so that the equation becomes

[ dfrac{P}{K−P}=C_1e^{rt}. label{eq30a}]

To solve this equation for (P(t)), first multiply both sides by (K−P) and collect the terms containing (P) on the left-hand side of the equation:

[egin{align*} P =C_1e^{rt}(K−P) [4pt] =C_1Ke^{rt}−C_1Pe^{rt} [4pt] P+C_1Pe^{rt} =C_1Ke^{rt}.end{align*}]

Next, factor (P) from the left-hand side and divide both sides by the other factor:

[egin{align*} P(1+C_1e^{rt}) =C_1Ke^{rt} [4pt] P(t) =dfrac{C_1Ke^{rt}}{1+C_1e^{rt}}. end{align*}]

The last step is to determine the value of (C_1.) The easiest way to do this is to substitute (t=0) and (P_0) in place of (P) in Equation and solve for (C_1):

[egin{align*} dfrac{P}{K−P} = C_1e^{rt} [4pt] dfrac{P_0}{K−P_0} =C_1e^{r(0)} [4pt] C_1 = dfrac{P_0}{K−P_0}. end{align*}]

Finally, substitute the expression for (C_1) into Equation ef{eq30a}:

[ P(t)=dfrac{C_1Ke^{rt}}{1+C_1e^{rt}}=dfrac{dfrac{P_0}{K−P_0}Ke^{rt}}{1+dfrac{P_0}{K−P_0}e^{rt}}]

Now multiply the numerator and denominator of the right-hand side by ((K−P_0)) and simplify:

[egin{align*} P(t) =dfrac{dfrac{P_0}{K−P_0}Ke^{rt}}{1+dfrac{P_0}{K−P_0}e^{rt}} [4pt] =dfrac{dfrac{P_0}{K−P_0}Ke^{rt}}{1+dfrac{P_0}{K−P_0}e^{rt}}⋅dfrac{K−P_0}{K−P_0} =dfrac{P_0Ke^{rt}}{(K−P_0)+P_0e^{rt}}. end{align*}]

We state this result as a theorem.

Solution of the Logistic Differential Equation

Consider the logistic differential equation subject to an initial population of (P_0) with carrying capacity (K) and growth rate (r). The solution to the corresponding initial-value problem is given by

[P(t)=dfrac{P_0Ke^{rt}}{(K−P_0)+P_0e^{rt}}].

Now that we have the solution to the initial-value problem, we can choose values for (P_0,r), and (K) and study the solution curve. For example, in Example we used the values (r=0.2311,K=1,072,764,) and an initial population of (900,000) deer. This leads to the solution

[egin{align*} P(t) =dfrac{P_0Ke^{rt}}{(K−P_0)+P_0e^{rt}}[4pt] =dfrac{900,000(1,072,764)e^{0.2311t}}{(1,072,764−900,000)+900,000e^{0.2311t}}[4pt] =dfrac{900,000(1,072,764)e^{0.2311t}}{172,764+900,000e^{0.2311t}}.end{align*}]

Dividing top and bottom by (900,000) gives

[ P(t)=dfrac{1,072,764e^{0.2311t}}{0.19196+e^{0.2311t}}.]

This is the same as the original solution. The graph of this solution is shown again in blue in Figure (PageIndex{6}), superimposed over the graph of the exponential growth model with initial population (900,000) and growth rate (0.2311) (appearing in green). The red dashed line represents the carrying capacity, and is a horizontal asymptote for the solution to the logistic equation.

Working under the assumption that the population grows according to the logistic differential equation, this graph predicts that approximately (20) years earlier ((1984)), the growth of the population was very close to exponential. The net growth rate at that time would have been around (23.1%) per year. As time goes on, the two graphs separate. This happens because the population increases, and the logistic differential equation states that the growth rate decreases as the population increases. At the time the population was measured ((2004)), it was close to carrying capacity, and the population was starting to level off.

The solution to the logistic differential equation has a point of inflection. To find this point, set the second derivative equal to zero:

[ egin{align*} P(t) =dfrac{P_0Ke^{rt}}{(K−P_0)+P_0e^{rt}} [4pt] P′(t) =dfrac{rP_0K(K−P0)e^{rt}}{((K−P_0)+P_0e^{rt})^2} [4pt] P''(t) =dfrac{r^2P_0K(K−P_0)^2e^{rt}−r^2P_0^2K(K−P_0)e^{2rt}}{((K−P_0)+P_0e^{rt})^3} [4pt] =dfrac{r^2P_0K(K−P_0)e^{rt}((K−P_0)−P_0e^{rt})}{((K−P_0)+P_0e^{rt})^3}. end{align*}]

Setting the numerator equal to zero,

[ r^2P_0K(K−P_0)e^{rt}((K−P_0)−P_0e^{rt})=0. onumber]

As long as (P_0≠K), the entire quantity before and including (e^{rt})is nonzero, so we can divide it out:

[ (K−P_0)−P_0e^{rt}=0. onumber]

Solving for (t),

[ P_0e^{rt}=K−P_0 onumber]

[ e^{rt}=dfrac{K−P_0}{P_0} onumber]

[ ln e^{rt}=ln dfrac{K−P_0}{P_0} onumber]

[ rt=ln dfrac{K−P_0}{P_0} onumber]

[ t=dfrac{1}{r}ln dfrac{K−P_0}{P_0}. onumber]

Notice that if (P_0>K), then this quantity is undefined, and the graph does not have a point of inflection. In the logistic graph, the point of inflection can be seen as the point where the graph changes from concave up to concave down. This is where the “leveling off” starts to occur, because the net growth rate becomes slower as the population starts to approach the carrying capacity.

Exercise (PageIndex{1})

A population of rabbits in a meadow is observed to be (200) rabbits at time (t=0). After a month, the rabbit population is observed to have increased by (4%). Using an initial population of (200) and a growth rate of (0.04), with a carrying capacity of (750) rabbits,

1. Write the logistic differential equation and initial condition for this model.
2. Draw a slope field for this logistic differential equation, and sketch the solution corresponding to an initial population of (200) rabbits.
3. Solve the initial-value problem for (P(t)).
4. Use the solution to predict the population after (1) year.
Hint

First determine the values of (r,K,) and (P_0). Then create the initial-value problem, draw the direction field, and solve the problem.

a. (dfrac{dP}{dt}=0.04(1−dfrac{P}{750}),P(0)=200)

b.

c. (P(t)=dfrac{3000e^{.04t}}{11+4e^{.04t}})

d. After (12) months, the population will be (P(12)≈278) rabbits.

Student Project: Logistic Equation with a Threshold Population

An improvement to the logistic model includes a threshold population. The threshold population is defined to be the minimum population that is necessary for the species to survive. We use the variable (T) to represent the threshold population. A differential equation that incorporates both the threshold population (T) and carrying capacity (K) is

[ dfrac{dP}{dt}=−rPleft(1−dfrac{P}{K} ight)left(1−dfrac{P}{T} ight)]

where (r) represents the growth rate, as before.

1. The threshold population is useful to biologists and can be utilized to determine whether a given species should be placed on the endangered list. A group of Australian researchers say they have determined the threshold population for any species to survive: (5000) adults. (Catherine Clabby, “A Magic Number,” American Scientist 98(1): 24, doi:10.1511/2010.82.24. accessed April 9, 2015, www.americanscientist.org/iss...a-magic-number). Therefore we use (T=5000) as the threshold population in this project. Suppose that the environmental carrying capacity in Montana for elk is (25,000). Set up Equation using the carrying capacity of (25,000) and threshold population of (5000). Assume an annual net growth rate of 18%.
2. Draw the direction field for the differential equation from step (1), along with several solutions for different initial populations. What are the constant solutions of the differential equation? What do these solutions correspond to in the original population model (i.e., in a biological context)?
3. What is the limiting population for each initial population you chose in step (2)? (Hint: use the slope field to see what happens for various initial populations, i.e., look for the horizontal asymptotes of your solutions.)
4. This equation can be solved using the method of separation of variables. However, it is very difficult to get the solution as an explicit function of (t). Using an initial population of (18,000) elk, solve the initial-value problem and express the solution as an implicit function of t, or solve the general initial-value problem, finding a solution in terms of (r,K,T,) and (P_0).

## Key Concepts

• When studying population functions, different assumptions—such as exponential growth, logistic growth, or threshold population—lead to different rates of growth.
• The logistic differential equation incorporates the concept of a carrying capacity. This value is a limiting value on the population for any given environment.
• The logistic differential equation can be solved for any positive growth rate, initial population, and carrying capacity.

## Key Equations

• Logistic differential equation and initial-value problem

(dfrac{dP}{dt}=rP(1−dfrac{P}{K}),P(0)=P_0)

• Solution to the logistic differential equation/initial-value problem

(P(t)=dfrac{P_0Ke^{rt}}{(K−P_0)+P_0e^{rt}})

• Threshold population model

(dfrac{dP}{dt}=−rP(1−dfrac{P}{K})(1−dfrac{P}{T}))

## Glossary

carrying capacity
the maximum population of an organism that the environment can sustain indefinitely
growth rate
the constant (r>0) in the exponential growth function (P(t)=P_0e^{rt})
initial population
the population at time (t=0)
logistic differential equation
a differential equation that incorporates the carrying capacity (K) and growth rate rr into a population model
phase line
a visual representation of the behavior of solutions to an autonomous differential equation subject to various initial conditions
threshold population
the minimum population that is necessary for a species to survive

## 9.6: Predator-Prey Systems - Mathematics

Part 1: Background: Canadian Lynx and Snowshoe Hares

In the study of the dynamics of a single population, we typically take into consideration such factors as the “natural" growth rate and the "carrying capacity" of the environment. Mathematical ecology requires the study of populations that interact, thereby affecting each other's growth rates. In this module we study a very special case of such an interaction, in which there are exactly two species, one of which -- the predators -- eats the other -- the prey. Such pairs exist throughout nature:

• lions and gazelles,
• birds and insects,
• pandas and eucalyptus trees,
• Venus fly traps and flies.

To keep our model simple, we will make some assumptions that would be unrealistic in most of these predator-prey situations. Specifically, we will assume that

• the predator species is totally dependent on a single prey species as its only food supply,
• the prey species has an unlimited food supply, and
• there is no threat to the prey other than the specific predator.

Very few such "pure" predator-prey interactions have been observed in nature, but there is a classical set of data on a pair of interacting populations that come close: the Canadian lynx and snowshoe hare pelt-trading records of the Hudson Bay Company over almost a century. The following figure (adapted from Odum, Fundamentals of Ecology , Saunders, 1953) shows a plot of that data.

To a first approximation, there was apparently nothing keeping the hare population in check other than predation by lynx, and the lynx depended entirely on hares for food. To be sure, trapping for pelts removed large numbers of both species from the populations -- otherwise we would have no data -- but these numbers were quite small in comparison to the total populations, so trapping was not a significant factor in determining the size of either population. On the other hand, it is reasonable to assume that the success of trapping each species was roughly proportional to the numbers of that species in the wild at any given time. Thus, the Hudson Bay data give us a reasonable picture of predator-prey interaction over an extended period of time. The dominant feature of this picture is the oscillating behavior of both populations.

1. On average, what was the period of oscillation of the lynx population?
2. On average, what was the period of oscillation of the hare population?
3. On average, do the peaks of the predator population match or slightly precede or slightly lag those of the prey population? If they don't match, by how much do they differ? (Measure the difference, if any, as a fraction of the average period.)

To be candid, things are never as simple in nature as we would like to assume in our models. In areas of Canada where lynx died out completely, there is evidence that the snowshoe hare population continued to oscillate -- which suggests that lynx were not the only effective predator for hares. However, we will ignore that in our subsequent development.

Part 2: The Lotka-Volterra Model

Vito Volterra (1860-1940) was a famous Italian mathematician who retired from a distinguished career in pure mathematics in the early 1920s. His son-in-law, Humberto D'Ancona, was a biologist who studied the populations of various species of fish in the Adriatic Sea. In 1926 D'Ancona completed a statistical study of the numbers of each species sold on the fish markets of three ports: Fiume, Trieste, and Venice. The percentages of predator species (sharks, skates, rays, etc.) in the Fiume catch are shown in the following table:

## What it means when the system contains higher-degree terms

Sometimes one or both equations will contain higher-degree terms. In the system

. x. is the predator population because . cxy. is positive, and . y. is the prey population because . -gxy. is negative. Because there are no higher-degree terms in the equation for . dy/dt. the prey population is only effected by the predators. On the other hand, because the equation for . dx/dt. contains the higher-degree term . bx^2. it means that the predator population is effected by the prey population, as well as another factor, like carrying capacity.

## Glossary

Some situations require more than one differential equation to model a particular situation. We might use a system of differential equations to model two interacting species, say where one species preys on the other. For example, we can model how the population of Canadian lynx (lynx Canadians) interacts with a the population of snowshoe hare (lepus americanis) (see https://www.youtube.com/watch?v=ZWucOrSOdCs).

A good historical data are known for the populations of the lynx and snowshoe hare from the Hudson Bay Company. This large Canadian retail company, which owns and operates a large number of retail stores in North America and Europe, including Galeria Kaufhof, Hudson's Bay, Home Outfitters, Lord & Taylor, and Saks Fifth Avenue, was originally founded in 1670 as a fur trading company. The Hudson Bay Company kept accurate records on the number of lynx pelts that were bought from trappers from 1821 to 1940. The company noticed that the number of pelts varied from year to year and that the number of lynx pelts reached a peak about every ten years. The ten year cycle for lynx can be best understood using a system of differential equations.

The primary prey for the Canadian lynx is the snowshoe hare. We will denote the population of hares by H(t) and the population of lynx by L(t), where t is the time measured in years. We will make the following assumptions for our predator-prey model.

If no lynx are present, we will assume that the hares reproduce at a rate proportional to their population and are not affected by overcrowding. That is, the hare population will grow exponentially,

The Lotka--Volterra system of equations is an example of a Kolmogorov model, which is a more general framework that can model the dynamics of ecological systems with predator-prey interactions, competition, disease, and mutualism. The equations which model the struggle for existence of two species (prey and predators) bear the name of two scientists: Lotka (1880--1949) and Volterra (1860--1940). They lived in different countries, had distinct professional and life trajectories, but they are linked together by their interest and results in mathematical modeling.

The predator–prey model was initially proposed by Alfred J. Lotka in the theory of autocatalytic chemical reactions in 1910. Lotka was born in Lemberg, Austria-Hungary, but his parents immigrated to the US. In 1925, he utilized the equations to analyze predator-prey interactions. Lotka published almost a hundred articles on various themes in chemistry, physics, epidemiology or biology, about half of them being devoted to population issues. He also wrote six books.

The same set of equations was published in 1926 by Vito Volterra, a mathematician and physicist, who had become interested in mathematical biology because of the impact by the marine biologist Umberto D'Ancona (1896--1964). Vito Volterra was born in Ancona, then part of the Papal States, into a very poor Jewish family. He attended the University of Pisa, where he became professor of rational mechanics in 1883. His most famous work was done on integral equations. In 1892, he became professor of mechanics at the University of Turin and then, in 1900, professor of mathematical physics at the University of Rome La Sapienza. His daughter Luisa married Umberto D’Ancona.

The predator-prey system of equations was later extended by many researchers, including C. S. Holling, Arditi--Ginzburg model, Rosenzweig-McArthur model, and some others.

The critical points of the Lotka--Volterra system of equations are the solutions of the algebraic equations

We may try to find the general solution of the Lotka--Volterra system of equations. From both equations, we get

Notice that the predator population, L, begins to grow and reaches a peak after the prey population, H reaches its peak. As the prey population declines, the predator population also declines. Once the predator population is smaller, the prey population has a chance to recover that the cycle begins again.

we can place two plots sude by side:

1. ( displaystyle 0 < a_1 delta_1 < K left( m_1 - delta_1 ight) ) and
2. ( displaystyle a_2 delta_2 > K left( m_2 - delta_2 ight) ) or ( b_2 le 1 . )

Mathematical analysis of the Beddington--DeAngelis system shows that there exist two equilibria (0,0) and ( left( frac , frac ight) = (4,1) , ) being globally stable in the interior of the first quadrant. The eigenvalues of the Jacobian matrix at the origin are ( lambda_1 =1 quadmboxquad lambda_2 =5 , ) and the eigenvalues of the Jacobian at the point (4,1) are given by ( - frac<1> <12>pm <f j>, frac> <12>. )

1. Barabas, Gyorgy, D'Andrea, Rafael, and Stump, Simon Maccracken, Chesson's coexistence theory, Ecological Monographs, 2018, https://doi.org/10.1002/ecm.1302
2. Batiha, K., Numerical Solutions of the Multispecies Predator-Prey Model by Variational Iteration Method, Journal of Computer Science, 2007, Vol. 3 (7): 523-527, 2007
3. Bayliss, A., Nepomnyashchy, A.A., Volpert, V.A., Mathematical modeling of cyclic population dynamics, Physica D: Nonlinear Phenomena, 2019, Volume 394, doi: 10.1016/j.physd.2019.01.010
4. Dellal, M., Lakrib, M., Sari, T., The operating diagram of a model of two competitors in a chemostat with an external inhibitor, Mathematical Biosciences, 302, No 8, 2018, 27--45.
5. Dimitrov, D.T. and Kojouharov, H.V., Complete mathematical analysis of predator–prey models with linear prey growth and Beddington–DeAngelis functional response, Applied Mathematics and Computation, 162, (2), 523--538, 2005.
6. Goel, N.S., Maitra, S.C., and Montroll, E., On the Volterra and other nonlinear models of interacting populations, Reviews of Modern Physics, 1971, Vol. 43, pp. 231-- https://doi.org/10.1103/RevModPhys.43.231
7. Hsu, S.B., Hubbell, S.P., and Paul Waltman, Competing predators, SIAM Journal on Applied Mathematics, 35, No 4, 1978, 617--625.
8. Hsu, S.B., Hubbell, S.P., and Paul Waltman, A mathematical theory for single-nutrient competition in continuous cultures of micro-organisms, SIAM Journal on Applied Mathematics, Vol 32, No 2, 1977, 366--383.
9. Hsu, S.B., Hubbell, S.P., and Paul Waltman, A contribution to the theory of competing predators, Ecological Monographs, Vol 48, No 3, 1978, 337--349.
10. May, R.M., Leonard, W.J., Nonlinear Aspects of Competition Between Three Species, SIAM Journal on Applied Mathematics, Vol. 29, Issue 2, pp. https://doi.org/10.1137/0129022
11. Molla, H., Rahman, M.S., Sarwardi, S., Dynamics of a Predator–Prey Model with Holling Type II Functional Response Incorporating a Prey Refuge Depending on Both the Species, International Journal of Nonlinear Sciences and Numerical Simulation, 2018, Vol. 20, No. 1, https://doi.org/10.1515/ijnsns-2017-0224
12. Olek. S. An accurate solution to the multispecies Lotka-Volterra equations, SIAM Review (Society for Industrial and Applied Mathematics), 1994, Vol. 36(3) (1994), 480–488 https://doi.org/10.1137/1036104
13. Ruan, J., Tan, Y., Zhang, C., A Modified Algorithm for Approximate Solutions of Lotka-Volterra systems, Procedia Engineering, 2011, Volume 15, 2011, Pages 1493-1497 https://doi.org/10.1016/j.proeng.2011.08.277
14. Scarpello, G.M. and Ritelli, D., A new method for the explicit integration of Lotka--Volterra equations, 2003, Divulgaciones Matematicas, Vol. 11, No. 1, pp. 1--17.
15. Seo, Gunog and Wolkowicz, Gail S. K., Sensitivity of the dynamics of the general Rosenzweig--MacArthur model to the mathematical form of the functional response: a bifurcation theory approach. Journal of Mathematical Biology, 76, No 7, 2018, 1873--1906.

## Preliminary

Letting P be prey density at time t, K the carrying capacity of the environment or segment buffering capacity, and assuming prey growth density over time obeys logistic pattern (P(1-P/K)) . If we consider predator handling rate is governed by Holling type II functional response of the form (Q/(a-P)) , where Q is the predator density at time t, we can write

where dP is the first order derivative, representing density variation of the prey in function of time, r is its intrinsic growth factor in the absence of the predator. This factor represent ratio of packets generated per time unit in respect to the ratio of prey packets present in the system.

(Q/(a+P)) is the functional response representing predator handling. a is a positive parameter denoting conversion rate of prey into resources.

Similarly, for the predator equation, we have

where dQ is the first order derivative representing predator density variation over time, β is its intrinsic growth factor in the absence of the prey. τ is the delay factor and tτ is the delaying terms denoting the fact that predator needs time to assimilate captured prey depending on system configuration. (P(t- au )/(a+P)) represents the functional response which corresponds to the assimilation efficiency. It is clear that predator growth over time is submitted to delay needed for processing predated preys and by the value of parameter a. A larger population of prey will ends up having a negative impact on predator population size, as the functional response will tend to one. Furthermore, in the absence of prey, as resources are limited and predator cannot grow exponentially, we assume the maximum population size the system can accommodate is K the carrying capacity (in terms of value to reach).

## 9.6: Predator-Prey Systems - Mathematics

Part 2.4: Equilibrium Points

Use the applet to experiment with trajectories for a predator-prey system, notice that there seems to be one set of initial conditions for which the trajectory consists of a single point. At such a point neither population function changes they are in equilibrium. For this reason such a point is said to be an equilibrium point .

Predator-Prey Direction Field and Trajectories

Return to your computer algebra worksheet. For our ongoing example, determine each of the coordinates of the equilibrium point of the system. Your calculations should be accurate to at least two significant digits.

(a) Consider the general predator-prey model.

If (x 0 ,y 0 ) is the equilibrium point, explain why x 0 and y 0 must satisfy

(b) Use the result in (a) to calculate the coordinates of the equilibrium point for the system with the coefficient values a =1 , b = 0.03 , c = 0.4 , and p = 0.01 . Compare the results of this calculation with your calculation in Step 1.

Let (x 0 ,y 0 ) be the equilibrium point for the general system

The lines x = x 0 and y = y 0 divide the first quadrant in the xy-plane into four subregions as indicated below.

For each subregion indicate whether x is increasing or decreasing and whether y is increasing or decreasing.

Predator-Prey Direction Field With Trajectories

For each of the four coefficients, a , b , c , and p , experiment with increasing and decreasing that coefficient while leaving the other three alone. When you have finished experimenting with a particular coefficient, return it to roughly the middle of its range. For each coefficient answer the following questions:

(a) How does the equilibrium point change as the coefficient changes?

(b) How does the pattern of trajectories change as the coefficient changes?

(c) Do the changes noted in (a) and (b) agree with your intuition about predator-prey interactions?

## Dynamics of a non-autonomous predator-prey system with Hassell-Varley-Holling Ⅱ function response and mutual interference

In this paper, we establish a non-autonomous Hassell-Varley-Holling type predator-prey system with mutual interference. We construct some sufficient conditions for the permanence, extinction and globally asymptotic stability of system by use of the comparison theorem and an appropriate Liapunov function. Then the sufficient and necessary conditions for a periodic solution of the system are obtained via coincidence degree theorem. Finally, the correctness of the previous conclusions are demonstrated by some numerical cases.

Citation: Luoyi Wu, Hang Zheng, Songchuan Zhang. Dynamics of a non-autonomous predator-prey system with Hassell-Varley-Holling Ⅱ function response and mutual interference[J]. AIMS Mathematics, 2021, 6(6): 6033-6049. doi: 10.3934/math.2021355

### Abstract

In this paper, we establish a non-autonomous Hassell-Varley-Holling type predator-prey system with mutual interference. We construct some sufficient conditions for the permanence, extinction and globally asymptotic stability of system by use of the comparison theorem and an appropriate Liapunov function. Then the sufficient and necessary conditions for a periodic solution of the system are obtained via coincidence degree theorem. Finally, the correctness of the previous conclusions are demonstrated by some numerical cases.

## Solution of Differential Equations

George Lindfield , John Penny , in Numerical Methods (Fourth Edition) , 2019

### 5.10.2 The Predator–Prey Problem

A system of differential equations which models the interaction of competing or predator–prey populations is based on the Volterra equations and may be written in the form

together with the initial conditions

The variables P and Q give the size of the prey and predator populations, respectively, at time t. These two populations interact and compete. K 1 , K 2 , C, and D are positive constants. K 1 relates to the rate of growth of the prey population P, and K 2 relates to the rate of decay of the predator population Q. It seems reasonable to assume that the number of encounters of predator and prey is proportional to P multiplied by Q and that a proportion C of these encounters will be fatal to members of the prey population. Thus the term CPQ gives a measure of the decrease in the prey population and the unrestricted growth in this population, which could occur assuming ample food, must be modified by the subtraction of this term. Similarly the decrease in the population of the predator must be modified by the addition of the term DPQ since the predator population gains food from its encounters with its prey and therefore more of the predators survive.

The solution of the differential equation depends on the specific values of the constants and will often result in nature in a stable cyclic variation of the populations. This is because as the predators continue to eat the prey, the prey population will fall and become insufficient to support the predator population which itself then falls. However, as the predator population falls, more of the prey survive and consequently the prey population will then increase. This in turn leads to an increase in the predator population since it has more food and the cycle begins again. This cycle maintains the predator and prey populations between certain upper and lower limits. The Volterra differential equations can be solved directly but this solution does not provide a simple relation between the size of the predator and prey populations therefore, numerical methods of solution should be applied. An interesting description of this problem is given by Simmons (1972) .

We now use Matlab to study the behavior of a system of equations of the form (5.28) applied to the interaction of the lynx and its prey, the hare. The choice of the constants K 1 , K 2 , C, and D is not a simple matter if we wish to obtain a stable situation where the populations of the predator and prey never die out completely but oscillate between upper and lower limits. The Matlab script below uses K 1 = 2 , K 2 = 10 , C = 0.001 , and D = 0.002 , and considers the interaction of a population of lynxes and hares where it is assumed that this interaction is the crucial feature in determining the size of the two populations. With an initial population of 5000 hares and 100 lynxes, the script e4s505.m uses these values to produce the graph in Fig. 5.13 .

Figure 5.13 . Variation in the population of lynxes (dashed line) and hares (solid line) against time, beginning with 5000 hares and 100 lynxes. Accuracy 0.005.

% x(1) and x(2) are hare and lynx populations.

simtime = input(ɾnter runtime ')

acc = input(ɾnter accuracy value ')

[t x] = ode23(fv,[0 simtime],initx,options)

xlabel('Time'), ylabel('Population of hares and lynxes')

For these parameters, there is a remarkably wide variation in the populations of hares and lynxes. The lynx population, although periodically small, still recovers following a recovery of the hare population.

## SIAM Journal on Applied Mathematics

A class of ODEs of generalized Gause type modeling predator-prey interaction is considered. The prey are assumed to exhibit a phenomenon called group defence, that is, predation is decreased or even eliminated due to the ability of the prey to defend or disguise themselves as their numbers increase.

Using the carrying capacity of the environment as the bifurcation parameter, it is shown that the model undergoes a sequence of bifurcations that includes a homoclinic bifurcation as well as a Hopf bifurcation. Conditions (that hold even in the case of no group defence) that ensure a subcritical Hopf bifurcation and also the spontaneous appearance of a semistable periodic orbit that splits into a pair (one stable and one unstable) of periodic orbits are given.

Ecological ramifications are considered. Unlike the classical model, sufficient enrichment of the environment combined with group defence leads to extinction of the predator (deterministically) for almost all initial conditions, providing strong support for the so-called “paradox of enrichment.”