19. The Lotka Volterra System
Section author: Malcolm Smith <msmith.malcolmsmith@gmail.com>
[status: first draft]
19.1. Motivation, Prerequisites, and Plan
Motivation
As discussed in the previous chapter, differential equations are a vital part of the way our world functions. They can be used to represent almost everything, from planetary motion to economics to electronics. In this chapter, we will focus on an ecological application of them, and how they can model the real world through abstraction.
Prerequisites
The 10-hour “Serious Programming” course
The “Simulating Differential Equations” mini course in Section 18
Plan
The plan for this mini course is to continue where the previous left off, with a focus on predator-prey relationships. If you missed the previous course, we will redo the introductory section, and then move on to the programming of the model.
19.2. A Brief Recap of Differential Equations
Before jumping in to the code, a description of what differential equations are and why they work can be found in Section 18.2.
19.3. The Lotka-Volterra System
One famous application of systems of differential equations is in predator-prey simulation, with varying degrees of complexity. The system we will look at is known as the Lotka-Volterra equations, and is a simple-but-interesting model. The model is based on the premise that there is some sort of closed ecosystem (an island, for example) with only two species: a predator and a prey. If these are wolves and deer, respectively, we can represent the Lotka-Volterra system with these equations:
These may look somewhat intimidating, but we can look at each term individually to understand what they mean. Before we begin, note that the greek letters a placeholders for some constants that we will assign later. The deer equation has two terms: a positive one proportional to the deer population, and a negative proportional to both the wolf and deer population. The first is simple to understand; if no other factors are present, the deer will grow exponentially. This has some unrealistic assumptions baked into it, but for now it’s a good enough approximation of reality. The second term is meant to model the rate of interactions between the wolves and deer. For the actual modeling, we will probably set \(\beta\) to 1, because it is a fair assumption that the deer will be killed in an interaction between the two.
For the wolf equation, the first term is also self-explanatory: without external influence, the wolf population will exponentially decay over time. This is because deer qualify as an external influence, so the wolves would starve without them. Then the last term represents the increase to the wolf population from gaining a single deer’s worth of meat. This is a hard thing to estimate, but it is nearly certain that the meat from one deer does not translate into more than one new wolf pup. Thus, it is likely \(\delta\) will be set below 1.
To simulate this, we can slightly modify our above code:
from matplotlib import pyplot as plt
import numpy as np
def main():
wolves = 1
deer = 1 ## initial conditions
wolf_pops = [wolves]
deer_pops = [deer]
step_count = 1500
step_size = 0.01
for _ in range(step_count):
wolves, deer = iterate(wolves, deer, step_size)
if wolves < 0.001:
wolves = 0
if deer < 0.001:
deer = 0
wolf_pops.append(wolves)
deer_pops.append(deer)
plt.plot(np.arange(step_count + 1) * step_size, deer_pops)
plt.plot(np.arange(step_count + 1) * step_size, wolf_pops)
plt.show()
def iterate(x, y, step_size):
"""step some parameters forward"""
alpha = 1.1
beta = -1
gamma = -0.4
delta = 0.5
x += (alpha * x + beta * x * y) * step_size
y += (gamma * y + delta * x * y) * step_size
return x, y
main()
Let us analyze what this code is doing. We have the initial conditions at the top, which can be chosen pretty much arbitrarily, as well as the length and granularity of the simulation. Then we update the populations the specified number of times, using the four constants described above. Finally, we plot these both against time, as can be seen below in Figure 19.3.1:
Figure 19.3.1 Two populations, oscillating regularly.
This is a fairly simple dynamic, where both populations seem to oscillate regularly. They are more complex than simple sine waves, though, as can be seen in the narrower peaks at the top of the orange line (wolf population).
Now that we have a simple model of this system, we can add more depth to it. There are several ways to do this, but the simplest is to simply add another population, representing vegetation. This is what the deer eat, and right now we are assuming an unlimited amount of it exists. This is an unrealistic assumption, so let us fix it:
import numpy as np
from matplotlib import pyplot as plt
def main():
grass = 0.1
deer = 0.8
wolves = 0.3 ## initial conditions
grass_pops = [grass]
deer_pops = [deer]
wolf_pops = [wolves]
step_count = 1500
step_size = 0.01
## g d w
rate_arr = [[-0.5, 0.3, 0], # g
[-1, 0.3, 0.7], # d
[0, -1, -0.5], # w
[3, -0.1, -0.5]] # base rate
for _ in range(step_count):
grass, deer, wolves = iterate([grass, deer, wolves], step_size, rate_arr)
if grass < 0.001:
grass = 0
if wolves < 0.001:
wolves = 0
if deer < 0.001:
deer = 0
grass_pops.append(grass)
deer_pops.append(deer)
wolf_pops.append(wolves)
plt.plot(np.arange(step_count + 1) * step_size, grass_pops)
plt.plot(np.arange(step_count + 1) * step_size, deer_pops)
plt.plot(np.arange(step_count + 1) * step_size, wolf_pops)
plt.show()
def iterate(pops, step_size, rate_arr):
"""step some parameters forward"""
rate_arr = np.array(rate_arr).T.copy()
d_pops = np.zeros_like(pops)
for i in range(len(pops)):
d_pops[i] += rate_arr[i][-1] * pops[i]
for j in range(len(rate_arr[i]) - 1):
d_pops[i] += rate_arr[i][j] * pops[i] * pops[j]
pops += d_pops * step_size
return pops
main()
Before discussing why this model is useful, we will look at what it does. Our previous model had 4 - or 2 squared - constants. This one has 12, which is more than 9 (3 squared). This will be explained in the next paragraph. For the most part, this program is otherwise the same. However, the iterate function now operates in a general manner, rather than the individual equations we had above. This means you could put in any number populations, provided you had the requisite array storing all the constants.
The reason this array is 3x4 instead of 3x3 is that we actually have four terms in each differential equation. Three of them are what you would expect; dependent on the current population, or the current population times a different one. The fourth is proportional to the square of the current population. The reason for this inclusion is that it help combat the inherently chaotic nature of systems of differential equations with more than two equations. This chaos is what’s at the root of many computationally difficult problems, such as weather forecasting or the three-body problem. Thus, the equations look like:
where N is the number of populations and \(P_i\) represents the ith population. The squared term here comes from when i equals j, and you have \(\beta _{ii} {P_i}^2\), where (as above) the greek letters represent a constant. In this case, the base rate row represents all the \(\alpha\) s, and the columns store all the constants that affect the particular population named at the top.
This model is more intersting, as there is more interaction between the three populations to be studied. We set the parameters so that the grass and the wolves’ interations have no bearing on the population, which is a fair assumption for this simple model. This model displays only limited periodicity, the reason for which will be explained below:
Figure 19.3.2 Three populations, oscillating wildly.
In this case, all three populations actually converge to some level, though it may appear as though they are oscillating randomly. This is because there are more terms in each differential equation, and thus more ways for all the terms to cancel each other out. In this case, an equilibrium is reached, with grass having the highest population, followed by deer, and then wolves, as can be seen below in Figure 19.3.3. This is a sign our model has somewhat realistic parameters, because the lower an animal is on the food chain, the more there are of it.
Figure 19.3.3 The same three populations that eventually stabilize.
19.4. Fractional Deer
So far, we have been assuming that any population can take any level it wants (provided that it is not negative). This allowed for smooth graphing, but does not make sense; What is a fractional deer (or wolf, or blade of grass)? In this case, we have been implicitly assuming that these numbers are in thousands, as evidenced by the less-than-zero check checking below 0.001. But beyond three digits of accuracy, these numbers wouldn’t make sense. So let us replace the less-than-zero checks:
if grass < 0.001:
grass = 0
if wolves < 0.001:
wolves = 0
if deer < 0.001:
deer = 0
with:
if grass < 0:
grass = 0
if wolves < 0:
wolves = 0
if deer < 0:
deer = 0
grass = round(grass, 3)
deer = round(deer, 3)
wolves = round(wolves, 3)
When you run this, it should look exactly like Figure 19.3.2, or very close. This is because, despite the unpredictability of these systems, this particular set of constants actually reaches a stable equilibrium, and small perturbations will not cause it to spiral into chaos. This may surprise you, because equilibriums in chaotic system tend to be unstable, like a ball balanced precariously at the top of a hill. However, this is why the Lotka-Volterra systems are a useful tool for modeling populations like these. They can model many types of stable patterns, similar to what might be found in nature.
The ability to not be phased by the fractional deer problem is a sort of litmus test for these models; if it gives significantly different results after a tiny tweak, it likely won’t work for real-world modeling. These particular sets of constants produced a good model, but playing around with other may produce results that do not have this quality.
For example, when the parameters are set like this:
## g d w
rate_arr = [[-0.5, -1, 0], # g
[0, -1, -2], # d
[-2.6, -1.6, -3], # w
[3, 4, 7.2]]# base rate
The program will exhibit a simple periodicity. However, the specific pattern will be significantly different depending on whether you round off or not, as seen below:
While they don’t seem too different, a closer look at the one on the right (the rounded off one) reveals a blockiness which is not found in the other graph. The difference is made much more stark by allowing the program to run longer (15000 steps):
Here is a significant difference between the two models. Interestingly, the model which levels out isn’t the base one, but the rounded one! For some reason, under this specific set of parameters, rounding to the thousands place curbs the oscillations to the point where they stabilize at one level. The unrounded model (which is actually rounded to 16 places because of the way python works) seems to encourage these oscillations, causing them to ramp up to absurd levels. Clearly, this system is inherently chaotic, and probably an unrealistic model of the real world (for one thing, there are a LOT more deer than grass).
There are other types of these equations (some that cause the extinction of a species, some that fluctuate randomly forever, etc), and they can be found in this paper.
19.5. Further Exploration
There are infinite combinations of possible values which could be the constant in a differential equation (or a system of them), and most of them either blow up to infinity or collapse to zero (for models like these, anyway). I would highly recommend taking a couple minutes to play around with those values, and see what interesting patterns emerge. Some values can actually create oscillation like we saw in the two population deer-and-wolf model, and it is a good thought experiment to try to think what real-world phenomenon could be represented by those particular parameters.
In addition, this is only one family of differential systems (we had families of solutions, and now families of systems… it’s families all the way down). There are thousands of other types, including those describing the quantum world, that are equally interesting to explore. Now that you understand what all (or most) of those symbols mean, there is not anything stopping you from understanding what all those equations represent. Given that these equations describe every effect in the universe we’ve ever observed, being able to read them is a very useful skill for pretty much everyone. In particular, a good family to practice on would be those describing Newtonian orbits, and the parallel ones describing the same orbits using the Hamiltonian (which is a topic worthy a chapter in its own right).
19.6. Conclusion
This has been an in-depth exploration of the Lotka-Volterra systems, which are a family of systems of differential equations that can be used to model the interactions between various populations in a closed system. We looked at the pros and cons of some approaches to modeling, and various patterns that can emerge from very simple starting rules. The goal of this chapter was both to help give an understanding of the specific system, and a general intuition about differential equations.